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1  Introduction 


Breast  cancer  is  the  most  frequently  diagnosed  malignancy  among  women  in  the  United 
States  [1].  In  1996  the  American  Cancer  Society  estimates  that  184,300  women  will  be 
newly  diagnosed  with  breast  cancer  and  that  44,300  will  die  from  the  disease  [1].  Breast 
cancer  accounts  for  31%  of  all  cancers  detected  and  17%  of  all  cancer  deaths,  and  ranks  as 
the  second  leading  cause  of  death  from  cancer  among  women  in  the  United  States  [1].  Five 
year  survival  rates  are  generally  very  high  (93%)  for  breast  cancer  staged  as  being 
localized,  falling  to  72%  for  regional  disease  and  only  18%  for  distant  disease  [2].  The  early 
detection  of  breast  cancer  is  clearly  a  key  ingredient  of  any  strategy  designed  to  reduce 
breast  cancer  mortality. 

Mammography’s  role  is  the  early  detection  of  breast  cancer.  Although  more  accurate  than 
any  other  modality,  existing  techniques  for  mammography  only  find  80  to  90%  of  the  breast 
cancers.  Moreover,  in  7  to  10%  of  cases,  the  cancer  will  not  be  visible  on  the  mammogram. 
It  has  been  suggested  that  mammograms  as  normally  viewed,  display  only  about  3%  of  the 
total  information  detected.  Perception  is  a  problem  particularly  for  patients  with  dense 
fibroglandular  patterns.  The  importance  of  diagnosis  of  breast  cancer  at  an  early  stage  is 
critical  to  patient  survival.  The  general  inability  to  detect  small  tumors  and  other  salient 
features  within  mammograms  motivates  our  investigation. 

The  goal  of  this  project  is  to  develop  an  interactive  diagnostic  tool  for  radiologists  that 
shall  refine  the  perception  of  mammographic  features  (including  lesions,  masses  and 
calcifications)  and  improve  the  accuracy  of  diagnosis.  By  improving  the  visualization  of 
breast  pathology  we  can  increase  the  chances  of  early  detection  of  breast  cancers  (improve 
quality)  while  requiring  less  time  to  evaluate  mammograms  for  most  patients  (lower  costs). 
We  are  investigating  a  methodology  for  accomplishing  mammographic  feature  analysis 
through  multiresolution  representations.  We  have  shown  that  efficient  (nonredundant) 
representations  may  be  identified  from  digital  mammography  and  used  to  reconstruct  and 
enhance  specific  mammographic  features  within  a  continuum  of  scale  space.  Such  focused 
reconstructions  may  complement  existing  modalities  and  allow  a  radiologist  to  examine 
interactively  diagnostic  features  within  a  selected  scale  space.  Similar  to  traditional  coarse 
to  fine  matching  strategies,  the  radiologist  may  first  choose  to  look  for  coarse  features  (e.g. 
dominant  mass)  within  low  frequency  levels  of  a  wavelet  transform  and  later  examine  finer 
features  (e.g.  microcalcifications)  at  higher  frequency  levels.  In  addition,  features  may  be 
extracted  by  applying  geometric  constraints  within  each  level  of  the  transform.  Choosing 
wavelets  (or  analyzing  functions)  that  are  simultaneously  localized  in  both  space  and 
frequency,  results  in  a  powerful  methodology  for  image  analysis. 
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A  major  reason  for  poor  visualization  of  small  malignant  masses  is  the  subtle  difference  in 
x-ray  attenuation  between  normal  glandular  tissues  and  malignant  disease  [3].  This  fact 
makes  the  detection  of  small  malignancies  problematical,  especially  in  younger  women  who 
have  denser  breast  tissue.  Although  calcifications  have  high  inherent  attenuation 
properties,  their  small  size  also  results  in  a  low  subject  contrast  [4].  As  a  result,  the 
visibility  of  small  tumors,  and  any  associated  microcalcifications,  will  always  be  a  problem 
in  mammography  as  it  is  currently  performed  using  analog  film. 

In  this  report,  we  describe  some  exciting  results  accomplished  during  the  past  year.  Below 
we  describe  in  executive  summary,  progress  related  to  the  specific  tasks  and  objectives 
stated  in  Phases  III  and  IV  of  our  original  Statement  of  Work. 

In  the  sections  below,  we  describe  in  overview,  our  processing  algorithms,  experimental 
methods  and  sample  results  obtained.  In  addition,  we  list  in  summary,  publications  of  our 
researchers  accomplished  during  the  past  year  of  our  investigation.  Finally,  we  include  our 
response  to  the  reviewers  comments  and  suggestions  regarding  contractual  and  technical 
issues  in  reference  to  last  year’s  report. 

1.1  Overview  of  Contents 

Generalizing  the  Discrete  Dyadic  Wavelet  Transform 

Traditional  orthogonal  and  biorthogonal  wavelet  representations  may  introduce  aliasing 
and  anisotropies  due  to  the  lack  shift  and  rotation  invariance.  In  screening  mammography 
intolerance  to  translation  and  rotation  is  clearly  not  desirable.  The  multiscale  spline 
derivative-based  transform  presented  in  this  section  is  a  redundant  representation  which 
does  not  introduce  such  artifacts. 

As  preliminary  results,  two  examples  of  image  fusion  are  presented.  Original  mammograms 
were  fused  with  their  globally  enhanced  counterparts  to  obtain  enhanced  images  having  an 
appearance  more  similar  to  the  original  mammograms.  The  fused  images  appears  easier  for 
radiologists  to  interpret  without  additional  training  on  wavelet  processing  techniques. 

A  Parallel  Algorithm  to  Support  Interactive  Wavelet  Processing  on  a 
Radiologist  Workstation 

We  have  designed  an  eflScient  parallel  algorithm  to  support  interactive  analysis  and 
multiscale  processing.  We  show  how  this  multiscale  block  algorithm  works  and  demonstrate 
the  amount  of  speedup  achieved  over  traditional  linear  convolution  schemes  and 
conventional  FFT  methods. 


12 


A  Continuous  Scale  Discrete  Wavelet  Transform  for  the  Detection  of  Spicular 
Masses  in  Mammograms 

We  introduce  a  simple  method  to  evaluate  continuous  wavelet  transform  (CWT).  We 
present  a  fast  decomposition  and  reconstruction  algorithm  to  perform  2-D  analysis  at 
arbitrary  scales.  Arbitrary  scale  analysis  is  shown  to  be  critical  for  optimal  feature 
detection  and  enhancement  in  mammograms  [5].  Our  algorithm  was  able  to  detect  a  mass 
that  could  not  be  seen  using  conventional  windowing  and  leveling  or  traditional  contrast 
enhancement  methods. 

Enhancement  of  Mammograms  from  Oriented  Information 

Mammograms  depict  all  the  significant  changes  in  breast  disease.  The  primary  radiographic 
signs  of  cancer  are  related  to  tumor  mass,  its  density,  size,  shape,  borders  and  calcification 
content.  These  features  may  be  extracted  according  to  their  coherence  and  orientation. 

Such  features  provide  important  visual  cues  for  radiologists  to  locate  suspicious  areas  easily. 
In  this  investigation,  an  artifact  free  enhancement  algorithm  based  on  multiscale  wavelet 
analysis  is  presented.  First  an  image  was  decomposed  using  a  fast  wavelet  transform 
algorithm.  At  the  same  time,  the  energy  and  phase  information  at  each  level  were 
determined  using  a  set  of  separable  steerable  filters.  Then  a  measure  of  coherence  within 
each  level  was  obtained  by  weighting  the  energy  with  the  ratio  of  projections  of  the  energy 
within  a  specified  window  onto  the  central  point  of  the  window  with  respect  to  the  total 
energy  within  each  window.  Finally,  a  nonlinear  operation,  integrating  coherence  and 
orientation  information,  was  applied  to  modify  the  transform  coefficients.  These  modified 
coefficients  were  then  reconstructed,  via  an  inverse  fast  wavelet  transform,  resulting  in  an 
improved  visualization  of  mammographic  features.  The  novelty  of  this  algorithm  lies  m  its 
detection  of  directional  features  and  removal  of  unwanted  perturbations.  Compared  to 
existing  multiscale  enhancement  approaches,  these  results  appear  more  familiar  to 
radiologists  and  naturally  close  to  the  original  mammogram. 

Multivoice  Undecimated  Wavelet  Transforms 

For  many  analysis/synthesis  applications  the  sampling  provided  by  the  traditional  wavelet 
transform  is  inadequate.  In  particular,  a  finer  grid  is  needed  when  computing  the  wavelet 
transform  for  applications  that  require  “good”  time-frequency  localization.  This  can  be 
achieved  by  computing  an  undecimated  wavelet  transform  and  the  addition  of  voices.  In 
this  section,  we  review  the  tools  that  will  allow  us  to  compute  an  exact  multivoice 
undecimated  wavelet  transform  for  an  arbitrary  wavelet.  Furthermore,  we  introduce  a 
wavelet  function  that  exhibits  nearly  optimum  time-frequency  localization  (the  sine-Gabor 
wavelet)  and  show  that  it  is  indeed  a  wavelet  frame.  However,  having  a  sine-Gabor  wavelet 
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with  nearly  optimum  time-frequency  localization  and  obtaining  a  tight  frame  at  the  same 
time  is  not  possible.  Both  constraints  are  important  because  they  allow  for  “good” 
time-frequency  localization  and  efficient  reconstruction,  respectively.  Obtaining  a  better 
(tighter)  frame  can  be  achieved  by  adding  voices.  This  not  only  eases  the  reconstruction 
but  is  also  a  requirement  for  applications  in  time-frequency  analysis.  Our  current  work  is 
concerned  with  extending  the  sine-Gabor  wavelet  to  multivoice  framework  and  with  the 
computation  of  multivoice  undecimated  wavelet  transforms.  Our  future  work  includes 
wavelet  reconstruction  in  the  multivoice/undecimated  framework  as  well  as  extending  our 
results  to  the  analysis  of  two  dimensional  data  for  mammography. 

Circular  Mass  Recognition  Based  on  the  Hough  Transform 

In  this  chapter  we  present  a  two-step  algorithm  for  the  recognition  of  circular  masses.  The 
first  step  uses  a  2D  Hough  Transform  for  the  detection  of  the  centers  of  possible  circular 
shapes  and  the  second  step  validates  their  existence  by  radius  histogramming.  The  2D 
Hough  Transform  described  in  this  section  makes  use  of  the  property  that  every  chord  of  a 
circle  passes  through  its  center.  We  present  sample  results  from  experiments  using  both 
synthetic  and  real  data  demonstrating  that  our  method  is  more  robust  to  noise  and 
complex  backgrounds  (typically  found  in  real  mammograms)  than  standard  gradient  based 
methods.  The  promise  of  the  method  is  demonstrated  with  its  application  on  a  digitized 
phantom  containing  four  circular  masses  and  on  a  digital  mammograms. 

1.2  Publications 

Below,  we  list  in  summary,  an  update  of  publications  accomplished  during  1996. 

1.  D.  Chen,  A.  Laine,  J.  G.  Harris,  and  W.  Huda,  “A  continuous  scale  discrete  wavelet 
transform  for  the  detection  of  spicular  masses  in  mammograms,”  summitted  to  IEEE 
Transactions  on  Signal  Processing,  February  1997. 

2.  D.  loannou,  W.  Huda,  A.  F.  Laine,  and  I.  Koren,  “Enhancement  of  computer  simulated 
masses  for  mammography  via  wavelet  analysis,”  summitted  to  IEEE  Transactions  on 
Medical  Imaging,  January  1997. 

3.  I.  Koren  and  A.  Laine,  “A  discrete  dyadic  wavelet  transform  for  multidimensional 
feature  analysis,”  in  M.  Akay  (Editor),  Time- Frequency  and  Wavelet  Transforms  in 
Biomedical  Engineering,  New  York,  NY:  IEEE  Press,  1997. 

4.  S.  Schuler  and  A.  Laine,  “Hexagonal  QMF  banks  and  wavelets,”  in  M.  Akay  (Editor), 
Time- Frequency  and  Wavelet  Transforms  in  Biomedical  Engineering,  New  York,  NY:  IEEE 
Press,  1997. 
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5.  J.  Fan  and  A.  Laine,  “Multiscale  contrast  enhancement  and  denoising  in  digital 
radiographs,”  in  A.  Aldroubi  and  M.  Unser  (Editors),  Wavelets  in  Medicine  and  Biology, 
Boca  Raton,  FL:  CRC  Press,  1996,  pp.  163-189. 

6.  A.  F.  Laine  and  X.  Zong,  “A  multiscale  sub-octave  wavelet  transform  for  de-noising  and 
enhancement,”  in  Wavelet  Applications  in  Signal  and  Image  Processing  IV,  Proceedings  of 
SPIE,  Denver,  CO,  August  6-9,  1996,  vol.  2825.  pp.  238-249. 

7.  A.  Laine,  W.  Huda,  D.  Chen,  and  J.  Harris,  “Segmentation  of  masses  using  continuous 
scale  representations,”  in  Digital  Mammography  ’96,  Proceedings  of  the  3rd  International 
Workshop  on  Digital  Mammography,  Chicago,  U.S.A.,  9-12  June  1996,  K.  Doi,  M.  L.  Giger, 
R.  M.  Nishikawa,  and  R.  A.  Schmidt,  Editors,  Amsterdam,  The  Netherlands:  Elsevier, 
1996,  pp.  447-450. 

8.  1.  Koren,  A.  Laine,  F.  Taylor,  and  M.  Lewis,  “Interactive  wavelet  processing  and 
techniques  applied  to  digital  mammography,”  in  Proceedings  of  the  IEEE  International 
Conference  on  Acoustics,  Speech  and  Signal  Processing,  Atlanta,  GA,  May  7-10,  1996,  vol. 
3,  pp.  1415-1418. 

9.  X.  Zong,  A.  F.  Laine,  E.  A.  Geiser,  and  D.  C.  Wilson,  “De-noising  and  contrast 
enhancement  via  wavelet  shrinkage  and  nonlinear  adaptive  gain,”  in  Wavelet  Applications 
III,  Proceedings  of  SPIE,  Orlando,  FL,  April  8-12,  1996,  vol.  2762,  pp.  566—574. 

10.  D.  Chen,  J.  G.  Harris,  and  A.  Laine,  “Automatic  scale  detection,”  in  Visual 
Communications  and  Image  Processing,  Proceedings  of  SPIE,  Orlando,  FL,  March  17-20, 
1996,  vol.  2727,  pp.  960-972. 

11.  G.  Qu,  W.  Huda,  and  C.  J.  Belden,  “Comparison  of  trained  and  untrained  observers 
using  subjective  and  objective  measures  of  imaging  performance,  in  Academic  Radiology, 
1996,  vol.  3,  pp.  31-35. 

1.3  Responses  to  Technical  and  Contractual  Issues 

We  believe  this  report  shows  significant  progress  with  respect  to  the  objectives  described  in 
Phase’s  III  and  IV  of  our  original  Statement  of  Work  (SOW).  We  acknowledge  the  reviews 
concerns  regarding  the  performance  of  our  methods  in  the  context  of  complex  backgrounds 
(typical  of  real  mammograms)  and  generating  false  positives.  Sections  1.1,  1.3,  1.4,  1.6 
describe  processing  techniques  that  are  artifact  free  and  significantly  minimize  the 
generation  of  false  positives  compared  to  previous  methods.  We  provide  experimental 
evidence  and  theoretical  support  to  justify  these  claims  in  each  case. 

In  addition,  with  regard  to  the  ROC  studies  mentioned  in  Phase  V  of  our  SOW,  we  include 
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below  a  description  of  our  ongoing  investigation.  It  is  expected  that  these  studies  shall  be 
complete  within  three  months  time. 

1.4  Study  in  Progress:  An  ROC  Comparison  Between  Digital 
Mammography  and  Screen-Film  Using  an  Anthropomorphic 
Breast  Phantom 

1.4.1  Introduction 

Several  studies  have  been  performed  comparing  the  performance  of  digital  mammography 
with  screen-film  mammography  with  ambiguous  results.  In  one  case,  the  conspicuity  of  low 
contrast  lesions  on  a  flat  grey  background  was  greater  using  digital  acquisition  and  display 
than  using  screen-film  combinations.  (1)  In  another  case,  a  similar  research  protocol  was 
used  with  simulated  masses  placed  on  an  anthropomorphic  breast  phantom,  resulting  in 
greater  conspicuity  using  the  screen-film  combination  than  with  the  digital  methods.  (2)  In 
spite  of  the  conflicting  results,  the  studies  provided  useful  for  the  following  reasons.  First, 
we  were  able  to  perfect  the  methodology  for  merging  simulated  masses  with  the  realistic 
background  of  the  anthropomorphic  breast  phantom,  and  second,  the  ROC  methodology 
developed  indicated  that  we  could  get  valid  results  using  individuals  with  a  wide  range  of 
experience.  The  abilities  to  generate  images  with  the  truth  established  and  to  evaluate  the 
images  without  the  requirement  for  sub-specialty  radiologist  readers  enable  us  to  use  a 
faster,  less  costly  approach  to  applying  established  ROC  methodology  to  the  evaluation  of 
our  processed  images. 

The  protocol  developed  by  Qu,  Huda,  et.al.  [6]  will  be  modified  slightly  to  be  used  to 
evaluate  steerable  wavelet  processing  of  digital  mammograms  as  described  in  Chapter  2.4. 
In  order  to  evaluate  and  exploit  the  power  of  the  steerable  wavelets,  two  shapes  of 
simulated  masses  will  be  used.  The  first  will  simulate  a  smooth,  mostly  round  mass,  while 
the  second  will  simulate  an  oblong  mass  with  rather  irregular  edges.  The  purpose  of  the 
study  is  to  evaluate  the  effectiveness  of  using  the  steerable  wavelets  on  different  shaped 
simulated  masses  in  the  context  of  complex  backgrounds.  Since  we  will  be  using  digital 
enhancement  methods,  only  digital  mammography  images  will  be  used  for  this  study. 

1.4.2  Method 
Phantoms 

Digital  images  will  be  generated  of  the  central  region  of  an  RMI  165  anthropomorphic 
phantom.  A  1  mm  thick  acrylic  base  will  be  positioned  on  top  of  the  RMI  165  phantom 
containing  a  3  x  3  array  of  1.67  cm  x  1.67  cm  regions  of  interest  (ROI).  Radiopaque 
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markers  will  be  located  in  the  four  corners  of  each  ROI  to  permit  these  regions  to  be 
identified  in  the  resultant  radiographic  images.  Each  image  will  be  5  cm  x  5  cm  in  size  and 
will  contain  9  identifiable  ROIs. 

The  test  objects  are  designed  to  simulate  masses  in  a  mammogram.  These  simulated 
masses  will  be  constructed  of  a  300  microns  (t)  film  polymer  and  positioned  in  the  central 
region  of  each  ROI  with  a  randomly  selected  orientation.  Two  types  of  simulated  masses 
will  be  made,  the  first  will  simulate  a  smooth,  mostly  round,  mass  with  a  diameter  of 
approximately  1.5  cm,  while  the  second  will  simulate  an  oblong  mass  with  rather  irregular 
edges  with  a  length  of  approximately  1.5  cm  and  a  width  of  approximately  0.5  cm.  Each  of 
the  masses  will  be  made  with  three  thicknesses,  5t,  4t,  and  3t  (six  total  mass  test  objects). 
Nine  radiographs  with  different  inserts  were  made  with  the  simulated  masses  located  in 
randomly  selected  ROIs  for  a  total  of  81  ROIs.  Approximately  half  of  the  ROIs  will 
contain  a  test  object  and  the  six  test  objects  will  be  evenly  distributed  randomly 
throughout  the  images.  In  addition,  the  test  objects  will  be  imaged  on  a  uniform 
background  to  be  used  to  demonstrate  the  expected  shapes  to  the  readers. 

Mammographic  Imaging 

Digital  mammographic  images  of  the  RMI  165  (+  insert)  will  be  generated  using  the 
LoRad  DSM  using  28  kVp/72  mAs  as  technique  factors.  The  images  will  be  transferred  by 
disk  to  a  sun  workstation  where  wavelet  processing  will  be  applied.  The  nine  images  (each 
with  nine  ROIs)  will  then  be  merged  to  form  a  single  image  with  a  156  pixel  black  border 
and  100  pixel  black  separators  between  the  nine  images  resulting  in  a  2048  pixel  x  2048 
pixel  image  (see  Figure  1).  Using  tools  developed  at  the  University  of  Florida,  the  image 
will  be  converted  to  a  DICOM  image  and  will  be  transferred  to  a  high  resolution 
workstation.  There  will  be  two  images  generated  with  phantom  data  to  contain 
unprocessed  image  data  and  wavelet  processed  image  data.  Each  512  pixel  block  will  either 
be  unprocessed  or  processed  and  the  blocks  will  be  randomly  placed  on  the  image.  The 
single  flat  background  image  with  the  six  test  objects  will  also  be  converted  to  a  2048  pixel 
X  2048  pixel  image  by  surrounding  the  image  with  a  768  pixel  black  border. 

Workstation  Display 

The  high  resolution  DICOM  workstation  developed  at  the  University  of  Florida  is  used  for 
diagnostic  interpretation  of  scanned  radiographs  and  computed  radiography.  The  megascan 
monitors  and  frame  buffers  are  driven  by  a  Sun  Sparc20  computer  with  image  display 
software  written  in-house.  Images  are  displayed  at  full  2048  x  2048  resolution  and  can  be 
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Figure  1:  Configuration  of  output. 


manipulated  with  zoom  /  roam,  window  /  level,  and  by  inverting  the  grey  scale. 

Reader  Performance 

Readers  will  be  shown  the  phantom  images  on  one  of  the  two  high  resolution  monitors  on 
the  megascan  with  the  flat  grey  image  demonstrating  the  test  objects  on  the  other.  By 
making  all  three  images  2048  pixels  x  2048  pixels,  we  are  assured  that  the  images  will  not 
be  subsampled  or  enlarged,  since  this  is  the  optimal  image  size  for  the  workstation.  We  will 
use  12  readers,  six  radiologists  (1  mammographer  and  5  senior  radiology  residents)  and  six 
scientists  experienced  in  reviewing  medical  physics  images  (two  basic  scientists  and  four 
graduate  students)  to  score  each  ROI.  The  readers  will  be  asked  to  review  each  ROI  and 
indicate  the  likelihood  of  the  ROI  containing  a  simulated  mass  using  a  scale  ranging  from  0 
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(test  object  definitely  absent)  to  10  (test  object  definitely  present).  A  rank  of  5  indicates 
that  there  is  a  50%  probability  of  the  simulated  mass  being  present. 

The  results  for  each  reader  and  processed  /  unprocessed  ROI  will  be  fitted  to  a  standard 
ROC  curve,  which  plots  the  True-Positive  Fraction  (TPF)  vs  False-Positive  Fraction  (FPF) 
as  the  detection  threshold  is  varied.  The  ROCFIT  software  package  (Metz,  University  of 
Chicago)  will  be  used  to  fit  the  acquired  True-Positive  Fraction  (TPF)  vs  False-Positive 
Fraction  (FPF)  data  and  to  obtain  the  resultant  area,  Az,  under  each  ROC  curve.  A 
paired  t-test  analysis  will  be  performed  to  assess  the  significance  of  the  differences  in  the 
area  under  the  ROC  curve. 
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Body 

2.1  Generalizing  the  Discrete  Dyadic  Wavelet  Transform 

2.1.1  Introduction 

Traditional  orthogonal  and  biorthogonal  wavelet  representations,  now  ubiquitous  in  image 
processing  and  computational  vision,  lack  shift  and  rotation  invariance,  and  may  introduce 
aliasing  and  anisotropies.  Such  artifacts  are  especially  problematic  for  medical  image 
analysis,  but  can  be  eliminated  by  the  use  of  redundant  representations. 

Since  a  mammogram  can  easily  be  translated  or  rotated,  it  is  not  desirable  that  the  resnlt 
of  processing  be  significantly  affected  by  positioning  during  acquisition  or  digitization.  The 
two-dimensional  transforms  presented  in  this  section  do  not  introduce  artifacts  due  to 
translation  and  rotation  invariance. 

First,  the  one-dimensional  discrete  dyadic  wavelet  transform  [7]  was  generalized  to 
higher-order  derivatives,  equipped  with  a  better  initialization  procedure,  and  used  for 
derivation  of  a  steerable  dyadic  wavelet  transform.  The  steerable  dyadic  wavelet  transform 
is  both  translation  and  rotation  invariant,  but  exact  reconstruction  is  problematic  via 
implementation  in  the  spatial  domain.  This  was  solved  by  a  spline-based  approximation 
which  was  used  for  a  multiscale  spline  derivative-based  transform. 

A  multiscale  spline  derivative-based  transform  was  implemented  in  the  spatial  domain  as  a 
filter  bank  consisting  of  x-y  separable  filters  only.  In  addition  to  not  introducing  artifacts, 
the  transform  enables  directional  analysis  of  images  across  dyadic  scales. 

As  a  sample  application,  two  examples  of  image  fusion  are  presented.  Original 
mammograms  were  fused  with  an  enhanced  version  to  obtain  an  enhancement  that  looked 
more  familiar  to  radiologists  than  mammograms,  enhanced  by  global  methods  previously 
developed. 

2.1.2  1-D  Discrete  Dyadic  Wavelet  Transform  Revisited 

In  our  previous  report,  we  presented  a  complete  description  of  discrete  dyadic  wavelet 
transform.  For  clarity,  we  summarize  below  definitions  of  the  continuous  and  discrete 
transform  [7]. 

The  dyadic  wavelet  transform  of  a  function  s(a;)  G  L‘^{R)  is  defined  as  a  sequence  of 
functions 

(1) 

where  Wras{x)  =  s  *  s(t)  dt,  and  i^rnix)  =  2-^^p{2-^x)  is  a  wavelet 

'0(x)  expanded  by  a  dilation  parameter  (or  scale)  2^.  To  ensure  coverage  of  the  frequency 
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axis  the  requirement  on  the  Fourier  transform  of  '0rn(^)  the  existence  of  Ai  >  0  and 
<  oo  such  that 

oc 

A^<  Y,  S 

m=—oo 

is  satisfied  almost  everywhere.  The  constraint  on  the  Fourier  transform  of  the  (nonunique) 
reconstructing  function  x(x)  is 

OO 

Y,  W2"a.)  x(2’"w)  =  1. 

m=—oo 

A  function  s(a;)  can  then  be  completely  reconstructed  from  its  dyadic  wavelet  transform 
using  the  identity 

OO 

=  X!  WmS*Xm{x), 
m=—oo 

where  Xm(^)  =  2— x(2-"^^). 

In  numerical  applications,  processing  is  performed  on  discrete  rather  than  continuous 
functions.  When  the  function  to  be  transformed  is  in  the  discrete  form,  the  scale  2’”  can  no 
longer  vary  over  all  tti  ^  Z .  Finite  sampling  rate  prohibits  the  scale  from  being  arbitrarily 
small,  while  computational  resources  restrict  the  use  of  an  arbitrarily  large  scale.  Let  the 
finest  scale  be  normalized  to  1  and  the  coarsest  scale  set  to  2^. 

The  smoothing  of  a  function  s(a;)  G  L^{R)  is  defined  as 

SjnS{x)  =  S*  (f>mix), 

where  (pmix)  =  with  me  Z,  and  (f>{x)  is  a  smoothing  function  (i.e.,  its 

integral  is  equal  to  1  and  (i>{x)  — >■  0  as  |a;|  — )•  oo). 

In  [7],  a  real  smoothing  function  (j){x)  was  selected,  whose  Fourier  transform  satisfied 

OO 

=  (2) 

m=l 

In  addition,  it  was  shown  that  any  discrete  function  of  finite  energy  s{n)  G  1^{Z)  can  be 
written  as  the  uniform  sampling  of  some  function  smoothed  at  scale  1,  i.e.,  s(n)  =  5'o/(n), 
where  f{x)  G  L'^iR)  is  not  unique.  Thus,  the  discrete  dyadic  wavelet  transform  of  s{n)  for 
any  coarse  scale  2^  was  defined  as  a  sequence  of  discrete  functions 


{SMf{n  +  s),  {Wmfin  +  s)}m6[l,M]}„eZ’ 


where  s  is  a  '0(x)  dependent  sampling  shift. 

The  above  initialization  s{n)  =  Sof{n)  is  rather  standard  in  the  discrete  wavelet  transform 
computation  [8],  although  it  yields  correct  results  (i.e.,  the  discrete  wavelet  transform  is 
equal  to  the  samples  of  its  continuous  counterpart)  only  when  s{n)  —  Sos{n).  Here,  we  will 
concentrate  on  wavelets  which  are  derivatives  of  spline  functions  and  this  will  lead  us  to  a 
simple  initialization  procedure  [9]  that  alleviates  the  above  problem. 

For  a  certain  choice  of  wavelets,  the  discrete  dyadic  wavelet  transform  can  be  implemented 
within  a  fast  hierarchical  digital  filtering  scheme.  Next,  we  shall  summarize  the  relations 
between  filters,  wavelets,  and  smoothing  functions. 

First,  let  us  introduce  a  real  smoothing  function  <f(x)  such  that  (2)  can  be  rewritten  as^ 

OO 

m~0 


and  let  us  set 


H^)  =  ^pix) 


(4) 


where  Pp{x)  denotes  a  central  B-spline  of  order  p.  With  the  choice  (4),  we  restrict  ourselves 
to  wavelets  which  are  spline  functions. 

Computing  (3)  for  the  finest  two  scales  shows  that 


ip^uj)  x(w)  =  I3p{u)  (p{oj)  -  Pp{2lv)  (p{2u:). 


(5) 


Pp{2im)  can  be  related  to  Pp{u))  by  expressing  Pp{2oj)  as  (cf.  Proposition  1  of  [9]) 


I3p{2uj) 


1  /  sin(a;) 

2P+1  I  sin  (f) 


p+1 


sin  (I) 


p-hl 


2 


and  using  E^=o 


;5p(2a;)  =  (cos  (1))^  (6) 

(Note  that  a  relation  similar  to  (6)  can  be  derived  for  integer  scales  provided  that  the 
dilation  parameter  and  the  order  p  are  not  both  even  [9].) 

Let  F(uj)  be  a  digital  filter  frequency  response  and  let  Fs(oj)  =  e^‘^^F(cj). 

^Note  that  the  sum  index  determines  the  range  of  scales  of  the  discrete  transform:  using  (2)  we  have 
and  x(2a;)  at  the  finest  scale  of  the  transform,  while  for  (3)  we  get  ipiuj)  and  x(w). 
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If  we  choose 


?^(a;)  =  G-s{^)^p{uj), 

(7) 

(p{2u)  =  Ls{lo)  (p{u)), 

(8) 

X{uj)  =  Ks{uj) 

(9) 

w)  =  e^‘""(^cos(|)y  , 

(10) 

where  s  G  {0,  is  a  filter  dependent  sampling  shift  needed  for  g{n),  l{n),  k{n),  and  h{n) 
to  be  FIR  filters,  and  insert  Equations  (6)-(10)  into  (5),  we  observe  the  relation  between 
frequency  responses  of  the  filters 


G{io)K{uj)  +  H{lo)L{u)  =  1. 


(11) 


Similar  to  orthogonal  and  biorthogonal  discrete  wavelet  transforms,  the  discrete  dyadic 
wavelet  transform  can  be  implemented  within  a  hierarchical  filtering  scheme.  To  derive 
such  a  digital  filtering  scheme,  let  us  assume  that  s{uj)  from  (1)  is  bandlimited  to  [-7r,7r]. 
Using  Shannon’s  sampling  theorem  [10]  and  (7)  in  the  definition  of  the  dyadic  wavelet 
transform  (1)  with  m  —  0  shows 

poo  ^ 

Wos{x)=  /  s{i)smc{t  -  i)  ^  g_s{m)Pp{x  -  t  -  m)  dt. 

^  i=—C)0  rn=—oo 

By  making  use  of  the  fact  that  the  cardinal  spline  functions  r]r{x)  tend  to  the  sine  function 
as  their  order  r  approaches  infinity  [11],  and  employing 

OO 

Vr{x)  =  b;:^i)l3r{x  -  i),  (12) 

i=—oo 

where  67  H«)  denotes  a  direct  B-spline  filter  [12],  we  can  write 


or,  by  using  (6)  and  (10), 

m— 1 

^{Wms{x)\^^^}  ~  5(0;)  jBp-|-r+i(^)  G^-s(2'"w)  ]^-ff-s(2*a;).  (13) 

Equation  (13)  entirely  specifies  the  discrete  dyadic  wavelet  transform  decomposition,  while 
the  reconstruction  follows  from  (5)— (10).  Three  levels  of  a  filter  bank  implementation  are 
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Figure  2:  Filter  bank  implementation  of  a  one-dimensional  discrete  dyadic  wavelet  transform 
decomposition  (left)  and  reconstruction  (right)  for  three  levels  of  analysis. 

shown  in  Figure  2.  (Note  that  the  initialization  is  the  same  as  the  one  proposed  in  [9]  and 
that  except  for  the  prefiltering  and  postfiltering,  this  scheme  is  implementing  “algorithme  a 
trous”  [13].)  Noninteger  shifts  at  scale  1  for  filters  with  s=|  are  rounded  to  the  nearest 
integer. 

We  will  now  perform  a  simple  experiment  which  will  illustrate  the  diiference  between  the 
implementation  of  the  discrete  dyadic  wavelet  transform  as  originally  proposed  in  [7]  (i.e., 
without  prefiltering  and  postfiltering)  and  the  one  from  Figure  2. 

Let  s(x)  —  sinc(x),  p  —  2,  and  ff(n)  —  25{n  H-  1)  —  25{n)  (this  particular  choice  for  p  and 
g{n)  results  in  the  same  wavelet  as  was  used  by  Mallat  and  collaborators  [14,  7]).  The 
dyadic  wavelet  transform  of  s{x)  at  a  scale  2™  (1)  in  the  frequency  domain  is  then 

W^s(uj)  =  G_,(2-a;)/?2(2'"a;)  rect(^)  .  (14) 

The  Fourier  transform  of  the  discrete  dyadic  wavelet  transform  of  s{n)  =  5(n)  at  a  scale  2™ 
using  spline  based  initialization  follows  from  (13) 

m— 1 

:F{WmS{n)}  =  Br+siio)  i7_,(2*n;),  (15) 

i~0 

while  the  one  using  the  algorithm  from  [7]  equals 

m~l 

B{Wms{n)}  =  G_,(2”^co)  J]  J7_,(2'a;).  (16) 

j=0 

In  Figure  3  a  comparison  of  magnitudes  of  (15)  and  (16)  versus  (14)  is  shown:  in  Figure 
3(a)  magnitudes  of  (14)  (solid)  and  (16)  (dashed)  are  plotted  for  m  €  {0, 1, 2, 3},  while  the 
dashed  curves  in  3(b)  represent  magnitudes  of  (15)  with  r  =  5. 

By  choosing  the  appropriate  order  r,  (15)  can  approximate  (14)  in  the  interval  [— tt,  tt] 
arbitrarily  good,  while  originally  proposed  (16)  has  troubles  at  finer  scales.  Mallat  and 
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Figure  3:  (a)  Fourier  transform  magnitudes  of  the  dyadic  wavelet  transform  of  s(a;)  =  sinc(a;) 
(solid)  and  the  originally  proposed  discrete  dyadic  wavelet  transform  [7]  of  s(n)  =  5{n) 
(dashed),  (b)  Fourier  transform  magnitudes  of  the  dyadic  wavelet  transform  of  s{x)  (solid) 
and  the  discrete  dyadic  wavelet  transform  using  quintic  splines  for  interpolation  of  s{n) 


Zhong  [7]  noticed  that  there  was  a  problem  with  their  discrete  transform  computation,  and 
introduced  a  set  of  constants  associated  with  the  discrete  transform  coefficients  at  dyadic 
scales.  They  chose  the  values  of  constants  such  that  the  transform  coefficient  modulus 
maxima  remained  constant  over  all  dyadic  scales  for  a  step  edge  input  signal.  In  relation  to 
Figure  3(a)  this  is  equivalent  to  multiplying  IF {Wmsin)}  by  a  distinct  constant  for  each  m. 
Clearly,  this  is  an  improvement  over  the  situation  depicted  by  Figure  3(a),  but  remains 
inferior  to  spline  based  initialization. 

Next,  we  choose  filters  in  the  filter  bank  implementation  of  the  discrete  dyadic  wavelet 
transform.  As  previously  mentioned,  we  are  interested  in  wavelets  which  are  derivatives  of 
spline  functions.  From  the  property  of  the  central  B-spline  functions  [15] 


dPp{.x) 

dx 


—  Pp-i 


it  follows  that  G{uj)  in  (7)  is  the  Fourier  transform  of  the  difference  operator  centered 
around  — s: 

G{lo)  =  (27  sin  ,  (17) 

where  d  is  the  order  of  the  derivative  and  the  sampling  shift  for  this  filter  is  s  = 

Since  H{u!)  was  already  given  by  (10),  the  remaining  two  filters  to  be  determined  are  L{(jj) 
and  K{uj).  Both  of  them  are  (as  is  true  for  (p{x)  and  x(^))  nonunique. 

If  we  choose  L(cu)  such  that  we  can  express  A(w)  in  terms  of  a  finite  geometric  series 
having  the  smallest  number  of  elements  for  an  arbitrary  p,  we  get 


and 


L{co)  =  ^  (-1)"^+^ 


(  L^J  ^ 

\  /  /CJX  \ 

(  COS  (  —  )  ) 

\  rn  , 

/  V  \2JJ 

K{u;) 


1 


dmod2 


(18) 


(19) 


where  [2;]  denotes  the  largest  integer  smaller  than  x,  the  sampling  shift  for  L{u)  is  the 
same  as  the  one  for  i7(a;)  (i.e.,  s  =  and  the  sampling  shift  for  K{uj)  is  the  same  as 

the  one  for  G{u)). 

Note  that  Equations  (10)  and  (17)-(19)  work  fine  from  the  mathematical  point  of  view, 
but  in  practice  the  reconstruction  may  become  cumbersome  when  both  p  and  d  are  large. 
(The  lengths  of  impulse  responses  h{n),  g{n),  l{n),  and  k{n)  are  p+2,  d+l, 
(p+l)(d-(d+l)mod2)+l,  and  pd+(p+l)(dmod2)+l,  respectively,  while  for  the  frequency 
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responses  of  the  decomposition  filters  we  observe  that  limp-^oo  and 

limd^oo(2j)"‘^  =  ^u{oJ  +  (2n+l)7r)  with  neZ.) 

It  is  also  worth  noting  that  K{lo)  is  a  lowpass  filter  when  p  is  even  (i.e.,  the  reconstruction 
function  x{^)  is  a  wavelet  only  for  p  odd). 

Tables  1,  2,  and  3  list  impulse  responses  of  the  four  filters  for  p  G  {0, 1, 2}  and  d  e  {1, 2, 3}, 
while  Figure  4  shows  wavelets  tp{x)  =  -  foi"  fh®  same  values  of  p  and  d.  Wavelets 

from  this  family  have  a  support  of  length  d+p+1,  regularity  order  p,  and  are  either 
symmetric  or  antisymmetric. 


Table  1:  Impulse  responses  h(n),  g{n),  l{n),  and  k{n)  for  p=0  and  d  G  {1,2,3}. 


n 

J.  . 

h{n) 

d-  1 

d-2 

K.  >  > 

d  =  3 

9{n) 

l{n) 

k{n) 

9{^) 

l{n) 

k{n) 

9{n) 

l{n) 

k[n) 

-2 

-1 

0 

1 

2 

0.5 

0.5 

1 

-1 

0.5 

0.5 

-0.25 

0.25 

1 

-2 

1 

0.5 

0.5 

-0.25 

1 

-3 

3 

-1 

-0.125 

0.625 

0.625 

-0.125 

0.0625 

-0.0625 

Table  2: 


Impulse  responses  h 


(n),  g{n),  Ijn),  and  k{n)  forp=l  and  d  G  {1,  2, 3} 


h(n) 

V'  ‘'7  5  ^  V' 

d 

=  1 

d 

=  2 

n 

l{n) 

9{n) 

k{n) 

9{n) 

k(n) 

-1 

0.25 

1 

-0.0625 

1 

-0.0625 

0 

0.5 

-1 

-0.3125 

-2 

-0.375 

1 

2 

0.25 

0.3125 

0.0625 

1 

-0.0625 

d  =  3 

n 

h{n) 

9{n) 

l{n) 

k{n) 

-3 

-0.015625 

-2 

1 

-0.09375 

0.00390625 

-1 

0.25 

-3 

0.265625 

0.04296875 

0 

0.5 

3 

0.6875 

0.1015625 

1 

0.25 

-1 

0.265625 

-0.1015625 

2 

-0.09375 

-0.04296875 

3 

-0.015625 

-0.00390625 

As  already  discussed,  wavelets  with  p  =  2  and  (i=l  from  a  family  of  wavelets  with  p  even 
and  d=l  were  used  in  [14,  7],  whereas  filters  with  p—l  and  d=2  from  a  family  of  filters 
with  p  odd  and  d=2  were  employed  in  [16,  17,  18].  Here  described  transform  puts  no 
restrictions  on  the  choice  of  p  or  d  whatsoever,  and  uses  a  better  initialization  procedure 
than  the  one  originally  proposed  in  [7]. 
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Table  3:  Impulse  responses  h{n),  g{n),  l{n),  and  k{n)  for  p=2  and  d  6  {1,2,3}. 

I  ^  ^  I  ^  I 


n 

h{n) 

9{n) 

l(n) 

k{n) 

9{n) 

l{n) 

k{n) 

-2 

-1 

0.125 

0.375 

1 

0.125 

-0.015625 

-0.109375 

1 

0.125 

-0.015625 

-0.125 

0 

0.375 

-1 

0.375 

-0.34375 

-2 

0.375 

-0.46875 

1 

0.125 

0.375 

0.34375 

1 

0.375 

-0.125 

2 

3 

0.125 

0.109375 

0.015625 

0.125 

-0.015625 

d  =  3 


n 

h{n) 

9{n) 

l{n) 

k{n) 

-4 

-0.001953125 

0.000244140625 

-3 

-0.017578125 

0.003662109375 

-2 

0.125 

1 

-0.0703125 

0.0263671875 

-1 

0.375 

-3 

0.085937 

0.0908203125 

0 

0.375 

3 

0.50390625 

0.13037109375 

1 

0.125 

-1 

0.50390625 

-0.13037109375 

2 

0.0859375 

-0.0908203125 

3 

-0.0703125 

-0.0263671875 

4 

-0.017578125 

-0.003662109375 

5 

-0.001953125 

-0.000244140625 

2.1.3  Steerable  Functions 

When  extending  the  transform  from  Section  2.1.2  to  two  dimensions,  one  of  the  first 
questions  that  come  to  mind  is  how  to  deal  with  the  fact  that  derivatives  can  be  defined  in 
any  direction  of  the  plane.  In  the  case  of  a  first  derivative  of  a  Gaussian,  one  can  simply 
compute  first  derivatives  of  a  Gaussian  in  x  and  y  directions  and  then  determine  the 
derivative  in  any  direction  from  these  two  directional  derivatives  [19].  For  higher  order 
derivatives  of  a  Gaussian,  Freeman  and  Adelson  [20]  showed  that  order+1  directional 
operators  are  needed  for  synthesizing  the  operator  at  any  orientation.  In  fact,  functions 
with  the  property  of  expressing  their  arbitrary  rotations  as  linear  combinations  of  some 
functions  are  not  limited  to  derivatives  of  a  Gaussian.  Below,  we  briefly  restate  some  of  the 
results  from  [20]. 

A  two-dimensional  function  is  called  “steerable”  if  its  rotations  generate  a  finite 
dimensional  space.  Rotations  of  a  steerable  function  /(r,  6)  can  therefore  be  expressed  as 

I 

/(r,S-9o)  =  E^W/i(>',»),  (20) 

i=l 

where  denotes  an  arbitrary  angle,  {ci(^o)}  stands  for  a  set  of  interpolating  functions. 
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{/j(r,  0)}  is  a  set  of  basis  functions,  and  r  =  \/x^  +y‘^  and  6  =  arg(x,  y)  are  polar  radius 
and  angle,  respectively. 

If  /(r,  9)  represents  a  filter  kernel,  the  result  of  filtering  with  a  rotated  filter  f{r,  6  —  ^o)  can 
be  computed  simply  by  {c,(^o)}  weighted  linear  combination  of  outputs  from  basis  filters 
{fi{r,  6)}.  Only  the  outputs  from  basis  filters  need  to  be  precomputed  and  then  the  output 
from  a  filter  rotated  by  any  angle  Oq  can  be  found  by  interpolating  between  them.  When  a 
large  number  of  rotations  of  a  template  filter  is  required,  the  above  scheme  can  lead  to 
substantial  savings  in  both  computational  cost  (time)  and  memory  requirements  (space). 
The  necessary  condition  for  a  function  /(r,  9)  to  be  steerable  is  that  /(r,  9)  is  bandlimited 
in  its  polar  angle; 

fir, 9)  =  anir)e^^f  (21) 

n=-N 


This  can  be  verified  by  inserting  (21)  into  (20)  and  by  assuming,  for  convenience,  that 
fi{r,9)  =  fir,  9-  9i): 

I 

a„(r)  c,(0o)a„(r)  (22) 

i=l 

where  {9i}  is  a  set  of  user  defined  angles  and  n  G  [-iV,  0].  ^  Since  only  nonzero  coefficients 
a„(r)  are  of  interest,  both  sides  of  (22)  can  be  divided  by  a„(r).  This  yields  a  matrix 
equation  from  which  interpolating  functions  Ci(0o)  can  be  determined: 


1  1  r  1  1  •••  1  cii9o) 

eJ^o  e.i02  ... 

...  ^jN9r  cii9o) 


(23) 


For  coefficients  =  0  the  rows  corresponding  to  each  n  were  removed  from  the  matrix 
formulation  shown  in  (23).  For  this  system  to  have  a  solution,  the  angles  {9i}  must  be 
chosen  such  that  the  columns  of  the  matrix  are  linearly  independent. 

In  [20]  they  proved  that  the  minimum  number  of  basis  functions  fi{r,9)  needed  to  steer 
fir,  9)  according  to  (20)  is  equal  to  the  number  of  nonzero  coefficients  a„(r)  in  the  Fourier 
series  expansion  (21). 

To  conclude  this  brief  description  of  steerability,  let  us  only  remark  that  functions  which 
are  not  steerable  (i.e.,  do  not  have  a  finite  number  of  terms  in  (21))  can  be  approximated 
with  steerable  functions  (a  singular  value  decomposition  was  employed  for  approximating 


^Note  that  the  constraints  are  the  same  for  n  G  [—N,  —1]  and  n  €  so  that  a  subset  of  all  possible 

values  for  n  G  [-N,  N]  can  be  taken. 
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such  functions  efficiently  [21]),  and  that  the  technique  of  expressing  transformed  versions  of 
a  function  as  linear  combinations  of  a  fixed  set  of  basis  functions  is  not  limited  to  rotations 
(extensions  to  translations  [22],  scalings  [21,  22],  and  general  transformations  [23]  have 
been  reported). 

2.1.4  Steerable  Dyadic  Wavelet  Transform 

Returning  to  the  question  from  the  beginning  of  Section  2.1.3,  the  answer  seems  obvious: 
one  needs  to  construct  a  steerable  analog  to  the  one-dimensional  transform  from  Section 
2.1.2.  Steerable  transforms  are  nothing  new — quite  a  few  [24,  25,  26,  22]  have  been  devised, 
some  of  them  [25,  26]  exactly  for  the  computation  of  directional  derivatives.  Here,  we  are 
not  interested  in  any  directional  derivatives:  we  want  to  use  derivatives  of  central  B-splines 
which,  as  the  order  of  B-splines  increases,  tend  to  derivatives  of  a  Gaussian. 

We  define  a  steerable  dyadic  wavelet  transform  of  a  function  s{x,y)  G  at  a  scale 

2”^,  m  e  Z,  as  [27] 

Wls{x,y)  =  s*i)i,{x,y),  (24) 

where  %l)l^{x^y)  denotes  'ipmiXiV)  rotated  by  6i,  tpm{^iy)  —  2  ^"*t^(2  ‘^x,2  ”*?/),  'ip{x,y)  is  a 
steerable  wavelet  that  can  be  steered  with  I  basis  functions,  and  di  =  with 
*  ^  2, . . .  ,  /}. 

Analogously  to  the  one-dimensional  case  (cf.  Section  2.1.2)  we  require  the  two-dimensional 
Fourier  plane  to  be  covered  by  the  dyadic  dilations  of  '0*(2”^Wj;,2  ^y)-  there  must  exist 
As  >  0  and  Bs  <  oc  such  that 

00  / 

A,  <  B,  (25) 

m=— oo  i=l 

is  satisfied  almost  everywhere. 

If  (nonunique)  reconstructing  functions  xlni^^v)  chosen  such  that  their  Fourier 
transforms  satisfy 

OO  I 

2'"a;,)  x'(2-a;„  2'"a;,)  =  1,  (26) 

m=— OO  i=l 

the  function  s{x,y)  may  be  reconstructed  from  its  steerable  dyadic  wavelet  transform  by 

OO  I 

E  (27) 

m=— oo  i=l 

where  xlni^^v)  denotes  Xmi^iy)  rotated  by  9i  and  Xmi^^y)  —  2  x(2  )r,2  y). 
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To  derive  an  algorithm  for  the  fast  computation  of  the  transform,  we,  similar  to  (3), 
introduce  two  smoothing  functions  such  that 

OO  I 

^{uJ^,U}y)ip{u^,LOy)  =  x"{2^UJa:,2^LOy).  (28) 

m=0  i=l 

We  choose  wavelets  that  are  steerable  analogs  to  the  one-dimensional  wavelets  from  Section 
2.1. 2;3 

/sinf^lX^’’’*^'''^ 

tp{Ur,  UJg)  =  {juJrCOs{LJe)Y  y  ^  J  ,  (29) 

where  cjr  =  +  ‘^y  ~  arg(<^x7^y)-  These  wavelets  can  be  steered  with  d-f-l  basis 

functions  (i.e.,  I  in  (20)  is  equal  to  d+1). 

Choosing 

^{2LOr)  =  HstiuJr)  (30) 

Y^{uJr,  iOe)  =  Gst{(^r,^^e  -  Oi)  4>{<^r),  (31) 

(p{2tjjY  ^ ^(a^r)j  (32) 

and 

X*(a;r,c^e)  =  Ksti<^r,<^e  -  Oi)  (p{iVr),  (33) 

and  using  (30)  through  (33)  with  (28)  computed  for  the  finest  two  scales,  we  obtain 

I 

^  ]  Gstif^T)  ^4>  Oi)Ksi{u}j-,U}(j,  6i)  Hst{j^r)Lst{u)r)  —  1.  (34) 

Z  =  1 

Setting  =  Pp{!^r),  and  employing  (6)  and  (29)  with  (30)  and  (31),  we  find  that 

HsM)  =  (35) 

and 

Gst{.oJr,tO0)  =  {cos{uje)YG-sYJr),  (36) 

where  H{(jo)  and  G{uj)  are  specified  by  (10)  and  (17),  respectively. 

By  inserting  (35)  and  (36)  into  (34),  the  missing  two  filters  can  be  chosen  as 

Lsti^r)  —  Ls{iOr)  (37) 


and 

Kst{oJr,<^e)  =  -^icos{u}e)YKs{LOr),  (38) 

where  L(a;)  and  K{oj)  are  given  by  (18)  and  (19),  respectively,  and  Cd  = 

®This  choice  can  be  viewed  as  an  extension  of  the  transform  presented  in  [28,  27,  29]. 
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2.1.5  Multiscale  Spline  Derivative-Based  Transform 

Let  us  pause  here  for  a  brief  assessment  of  the  two-dimensional  steerable  transform  derived 
so  far.  We  have  chosen  steerable  wavelets  (29)  which  are  equal  to  d-th  order  derivatives  of 
circularly  symmetric  spline  functions  in  the  direction  of  x-axis  (note  that  knots  for  these 
splines  are  circles)  and  we  have  laid  a  foundation  for  filter  bank  implementations  in  (34). 
An  obvious  shortcoming  of  this  scheme  is  the  fact  that  none  of  the  filter  kernels  obtained 
from  (35)  through  (38)  is  compactly  supported  on  the  rectangular  grid.  For 
implementations  using  digital  filters,  one  is  therefore  forced  to  approximate  these  frequency 
responses,  and  by  doing  so,  (34)  may  not  hold  anymore.  Filters  in  filter  bank 
implementations  of  steerable  pyramids  described  in  [25,  26,  22],  for  example,  were  designed 
by  using  various  techniques  for  approximating  desired  frequency  responses.  None  of  the 
reported  filter  banks  achieved  perfect  reconstruction  and  all  filter  kernels  were 
nonseparable.  Here,  we  will  take  a  different  approach.  We  will  approximate  the  wavelets  in 
(29)  in  a  way  that  will  lead  to  x-y  separable  filters  in  a  perfect  reconstruction  filter  bank 
implementation  of  the  transform  such  that  the  quality  of  approximation  will  improve  by 
increasing  the  order  of  B-splines. 

Let  us  approximate  wavelets  from  (29)  with 

=  (39) 

Based  on  the  fact  that  B-splines  tend  to  a  Gaussian  as  their  order  increases,  it  is  easy  to 
see  that  both  wavelets  (29)  and  (39)  converge  to  the  same  functions  (i.e.,  d-th  order 
derivatives  of  the  normalized  Gaussian  in  the  direction  of  a;-axis)  as  p  -)■  oo. 

In  order  to  steer  wavelets  'ip{x,y)  given  by  (39)  (note  that  steering  will  be  only 
approximate,  since  these  wavelets  are  not  steerable),  we  need  to  find  basis  functions  that 
will  approximately  steer  ip{x,y).  Computing  rotations,  as  we  did  in  (24),  is  not  practical 
here,  because  arbitrary  rotations  of  (39)  cannot  be  expressed  exactly  in  terms  of  central 
B-spline  functions.  Instead,  we  take  advantage  of  the  property  of  circularly  symmetric 
functions  that  rotations  of  their  directional  derivatives  are  equal  to  directional  derivatives 
in  rotated  directions; 

^  fd'^Qcix,y)\  d^Qo{x,y) 

dn^  J”  dfil  ’ 

where  Ue,  stands  for  rotation  by  do,  =  nVgc{x,y),  Qc{x,y)  is  a  circularly  symmetric 

function,  and  ngp  denotes  vector  n  =  (cos  d,  sin  d)  rotated  by  do. 

Let  us  choose 

Q{x,y)  =  f3p+d{x)Pp+diy), 
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which  is  approximately  circularly  symmetric  function  for  higher  order  splines.  A  rotation  of 
tp{x,y)  =  from  (39)  by  6o  can  therefore  be  approximated  by 


d^gjx,  y) 
dn^ 


d 

E 

2=0 


d 


id^  ^Pp+d{x)  d^l3p+d{y) 


dx^~ 


dy^ 


(40) 


where  ft  =  (cos^o^sin^o)  =  (nx.ny). 

Note  that  in  case  of  Gaussian,  which  is  both  x-y  separable  and  circularly  symmetric,  (40) 
becomes  exact  (e.g.,  for  Q{x,y)  =  9o  =  -0,  and  d  =  {2,4},  we  obtain,  up  to  a 

scaling  factor,  x-y  separable  basis  set  for  the  second  and  fourth  derivative  of  Gaussian  from 
Tables  III  and  VII  of  [20]). 

Having  found  a  set  of  basis  functions  (40)  that  approximately  steer  wavelets  (39),  we  want 
to  construct  a  transform  such  that  Equations  (24)  through  (28)  will  be  valid  (superscript  i 
must  be  viewed  now  as  an  index,  rather  than  rotation  by  9i).  In  frequency  domain,  we  can 
express  basis  functions  from  (40)  as 


=  GiJ{uJx)Gl,{(Jy)pp+iM^p+d-iM,  i  =  0, 1, . . .  ,d. 


(41) 


where  G^{u)  is  given  by  (17)  and  G^{(jj)  —  1. 

Choosing  appropriate  'jd'(ujx,ujy)  to  obtain  a  relation  needed  for  the  filter  bank 
implementation  of  the  transform  is  more  complicated  than  in  one  dimension.  Since  we 
could  not  find  a  general  solution  for  arbitrary  d,  we  solve  each  case  separately.  Below,  we 
present  solutions  for  the  first  four  orders.  When  d  <  2,  we  impose 

(p{uix,<^y)  =  ^{cOx,(^y)  =  ^pi^x) Ppi^y) ,  &  Constraint  analogous  to  the  one  from  Section  2.1.2. 
For  d  =  1,  we  write  similar  to  [7] 

x\LOx,t0y)  =  Kl{(Vx)Ti{Uy)$p{uJx)^p-l{u^y),  (42) 

X^{u}x,UJy)  =  Ti{uJx)Kl{oJy)Pp^i{uJx)Pp{<^y),  (43) 


where 

T,(w)=i(n-iffMn 

and  K^{uj)  is  given  by  (19). 

Computing  (28)  for  the  finest  two  scales  and  inserting  (6),  (41),  (42),  and  (43)  yields  a 
relation  between  frequency  responses 

G^{uJx)K^{‘^x)Ti.{^y)  +  Ti{LL>x)G^{u)y)K^{u!y)  +  \H{u>x)H{<^y)\‘^  =  1- 
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For  d  =  2,  we  choose 


(44) 

(45) 

(46) 


—  Kg{uJx)T2{^y)Pp{^x)Pp~2{^y)^ 

X^i^xj^y)  ~  Kl{uJx)Kl{Uy)^p-i{cOx)^p-l{^y)^ 

X^(u)x^UJy)  =  T2{u)x)Kg{uJy)^p^2{^x)Pp{^y)i 

where 

T2{uj)  =  \H{uj)\^,  (47) 

Using  (6),  (41),  and  (44)  through  (46)  with  (28)  results  in 

G‘^{u)x)K'^{(^x)T2{^y)  G^{^x)K^{^x)G^{^y)^  (^y)  +  72(^2;)^  i^y)^  (^?/)"^ 

+|i/(a-Jff(w,)r  =  1. 

For  orders  d  >  2,  we  require  i(u„u,)  =  /3p(wi)4(w,)  and  ip{u„ui,)  =  where 


(P(lo)  is  specified  by  (8)  and  (18). 

For  0?  =  3,  we  choose  reconstructing  functions 

-^{(jJx^OJy)  =  —Kl((jL}x)Kl{u)y)Vz{oJx)yz{^y)^-l{.^x)^-2{^y)-,  (^9) 

X^(a;a;,<^2,)  =  -Kl{iOx)Kl{uJy)Vz{iOx)Vz{iOy)<P-2{oJx)<P-l{f^y)-,  (59) 

=  Kl{uJy)^-z{u)x)<PM,  (51) 

where 

V^(uj)  =  i(l  -  |//(n.)p).  (62) 


and  €  L^{R)  denotes  a  function  such  that  (p{uj)  =  Pi-i{u>)(p-i{io),  i  E  N. 

Employing  (41),  (6),  (8),  and  (48)  through  (51)  with  (28)  yields  a  relation 

G^{uJx)K^{tOx)  -  G‘^{uJx)K^{i^x)Vz{Ux)G^{oJy)K^i^^y)V3{Uy)~ 
-G\cOx)K\ux)Vz{ux)G\uy)K\ujy)Vz{u;y)  +  G^Uy)K^{u;y)+ 

-\-H(u)x)L(jjJx)H[(Jy')L[LOy)  1, 

where  L{u})  is  specified  by  (18). 

For  d  =  4,  our  choices  are 


x\Ux,C0y)  =  Kt{i0x)T2{uJy)(p{uJx)<P-AM,  (53) 

X^{uJx,UJy)  =  K^sMKl{uJy)(p-i{uJx)^-z{(^y),  (54) 

x\uJx,COy)  =  -K^A^x)K^sMM^x)V4{uJy)^-2M0-2M,  (55) 

X^{cJx,UJy)  =  Kl{uJx)K^{(^y)(P-3{oJx)^-i{‘^y),  (59) 

X^(L0x,^^y)  =  T2{uJx)K^i(^y)^-i{^x)0{^^y),  (5"^) 
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where 


Vi{u,)  =  l-\H{Lo)\\  (58) 

Using  the  above  functions  (53)  through  (57),  (41),  (6),  and  (8)  in  (28)  computed  for  the 
finest  two  scales  gives 


G^{,^x)K'^{!-^x)T2{u}y)  +  G'^  (ux)K^  {uJx)G^  {u)y)K^  {ojy)— 

—  G  +  G^  {lOx)K^  {fjJx)G^  {u)y)K^  {ujy)  + 

+T2{uJx)G^{Uly)K‘^{u}y)  +  H{(jJx)L{Ux)H{(jJy)L(u)y)  =  1. 


Here,  we  have  even  more  freedom  for  choosing  the  reconstructing  functions  than  in  one 
dimension.  The  above  functions  for  d  =  {2,3,4}  were  found  by  trying  to  reproduce  the 
one-dimensional  transform  from  Section  2.1.2  as  much  as  possible.  All  decomposition  filters 
G^(a;)  were  first  paired  with  corresponding  reconstruction  filters  K^{oo),  and  then  all  other 
potential  digital  filters  were  specified  as  polynomials  in  We  inserted  thus  specified 

filters  into  a  relation  between  their  frequency  responses  and  solved  for  the  unknown 
polynomial  coefficients.  Since  we  allowed  more  filters  with  higher-degree  polynomials  than 
expected  in  the  solution,  the  resulting  system  of  linear  equations  was  under  determined. 
This  allowed  enough  freedom  for  removal  of  undesired  digital  filters  and  for  balance 
between  degrees  of  polynomials. 

The  described  procedure  for  determination  of  reconstructing  functions  and  filters  involves 
quite  a  lot  of  heuristics  to  obtain  the  appropriate  solution  from  the  underdetermined  linear 
system.  Unfortunately,  we  are  not  aware  of  any  systematic  way  (aside  from  numerical 
optimization,  which  may  be  pretty  cumbersome)  to  obtain  solutions  comparable  to  the 
ones  above. 

Next,  we  will  derive  a  filter  bank  implementation  of  the  transform.  Similar  to  the 
one-dimensional  case  from  Section  2.1.2,  we  assume  a  bandlimited  input  signal: 
s{u>x,ojy)  =  0  for  \u}x\  >  TT  or  \ujy\  >  tt.  Using  Shannon’s  sampling  theorem  in  two 
dimensions  [30]  with  (24)  and  basis  functions  from  (41),  we  can  write 


nOQ  OO  OO 

s{ix,  iy)  sinc(ta;  -  4)  sinc(4  -  iy)- 

ix=—ooiy=—oo 


OO  OO 

■  9t'J'{rnx)l5p+i{x  -tx-  rux)  g''_^{my)(5p+d^i{y  -  ty  -  my)  dtx  dty, 

mx=—oo  rny=—oo 


where  z  =  0, 1, . . .  ,  d  as  in  (41). 
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Again,  we  approximate  sine  functions  with  r-order  cardinal  splines,  then  use  (12),  and  get 

~  S{ux,ujy)  B;\uy)  Bp+r+i+iM- 

m— 1 

■Bp+,+j_j+iK)  Gl7(2”a-,)  G*_.(2'“a-,)  J]  (59) 

71=0 

Using  (59)  with  an  approximation  Bp^r+i+i{^)  —  Bp^r{i^)Bi{uj),  we  can  obtain  a  filter 
bank  implementation  of  the  transform  decomposition.  The  reconstruction  part  follows  from 
(28),  (41),  and  reconstructing  functions  for  distinct  values  of  d.  Figure  5  shows  filter  bank 
implementations  of  a  multiscale  spline  derivative-based  transform  for  d  =  {1, 2, 3, 4}.  For 
d  =  1,  we  recognize  (except  for  the  prefiltering  and  postfiltering)  the  filter  bank 
implementation  of  a  two-dimensional  discrete  dyadic  wavelet  transform  from  [7].  For  d  —  2, 
however,  our  transform  differs  from  the  filter  bank  presented  in  [17];  second  derivative  is 
computed  only  in  the  directions  of  x  and  y-axis  in  [31,  17],  which  is  not  enough  for 
steering.  Although  not  particularly  appropriate  for  orientation  analysis,  such  a  scheme  may 
still  efficiently  approximate  Laplacian  of  Gaussian  across  dyadic  scales. 

A  transform  similar  to  the  one  described  in  this  section,  was  presented  in  [25,  26,  22]. 

Their  filter  bank  implementation  is  less  redundant  (downsampling  is  used  on  the  lowpass 
branch,  while  simultaneously  keeping  aliasing  negligible  by  using  a  filter  with  insignificant 
amount  of  energy  for  \u>r\  >  |)  and  uses  reconstruction  filters  with  same  magnitude 
frequency  responses  as  the  decomposition  ones — a  possible  advantage  for  certain 
applications.  They  have,  on  the  other  hand,  little  control  over  the  function  from  which 
derivatives  are  computed  (to  obtain  a  d-th  derivative,  they  multiply  a  circularly  symmetric 
filter  by  {—j  cosOy  with  all  filters  being  obtained  by  recursive  minimization  of  a  weighted 
combination  of  constraints),  filter  bank  does  not  have  perfect  reconstruction,  and  none  of 
the  filters  is  x-y  separable. 

2.1.6  Image  Fusion  Application 

A  new  direction  in  our  visualization  methodology  (Phase  IV  statement  of  work)  involves 
the  use  of  image  fusion.  A  problem  with  our  early  algorithms  using  multiscale  global 
enhancement  was  that  they  produced  results  which  were  somewhat  unfamiliar  looking 
compared  to  traditional  methods  of  enhancement  and  the  original  mammogram.  Examples 
of  two  such  cases  are  shown  in  the  middle  of  Figures  6(a)  and  6(b):  textures  of  fibrous 
tissues  extracted  from  the  modulus  of  transform  coefficients  and  multiscale  histogram 
equalization  [32].  These  representations  did  an  outstanding  job  of  extracting  salient 
features  in  mammograms  but  were  difficult  for  radiologists  to  interpret  without  additional 
training. 


37 
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Figure  5:  Filter  bank  implementation  of  a  multiscale  spline  derivative-based  transform  for 
m  G  [0,M  -  1]:  (a)  Prefiltering,  (b)  postfiltering,  (c)  decomposition  and  (d)  reconstruction 
modules  for  d  =  1,  and  (e)  decomposition  and  (f)  reconstruction  modules  for  d=2. 
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(i)  (j) 


Figure  5:  Continued:  (g)  Decomposition  and  (h)  reconstruction  modules  for  d  =  3,  and 
(i)  decomposition  and  (j)  reconstruction  modules  for  d  =  4.  Decomposition  modules  are 
recursively  applied  at  the  locations  of  the  filled  circles. 
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By  integrating  the  original  mammogram  with  (one  or  more)  enhanced  versions  we  can 
produce  (synthesize)  an  image  which  is  more  familiar  (and  appealing)  to  mammographers 
and  general  radiologists.  Thus,  fusion  plays  an  important  role  in  avoiding  a  “paradigm 
shift”  in  current  protocols. 

Traditionally,  image  fusion  combines  different  aspects  of  information  from  the  same 
imaging  modality  or  from  distinct  imaging  modalities  [33]  and  can  be  used  to  improve  the 
reliability  of  a  particular  computational  vision  task  or  to  provide  a  human  observer  with  a 
deeper  insight  about  the  nature  of  observed  data. 

Our  fusion  method  is  based  upon  a  multiscale  spline  derivative-based  transform  [28]  which 
enables  multiscale  image  processing  along  arbitrary  orientations.  Steerable  dyadic  wavelet 
transforms  share  many  properties  with  multiscale  representations  used  in  image  fusion, 
such  as  pyramids  [34,  35]  and  traditional  wavelet  analysis  [36,  37].  However  an  important 
difference  is  translation  and  rotation  invariance,  and  the  absence  of  aliasing  in  its  filter 
bank  implementation.  Both  aliasing  and  translation  non-invariance  can  be  sources  of 
unwanted  artifacts  in  a  fused  image  [38]. 

For  images  to  be  fused,  filters  from  the  filter  bank  implementation  of  a  multiscale  spline 
derivative-based  transform  are  used  in  quadrature  pairs  (i.e.,  with  their  Hilbert  transform 
counterparts) .  Quadrature  pairs  of  filters  steered  to  some  arbitrary  angle  are  used  to 
determine  the  local  oriented  energy  defined  as  the  sum  of  the  squared  outputs  from  each 
filter  of  the  quadrature  pair.  Local  dominant  orientation  (i.e.,  the  angle  that  maximizes  the 
local  oriented  energy)  is  determined  at  each  level  and  position,  filters  are  steered  to  the 
local  dominant  orientation,  and  local  oriented  energies  are  compared.  The  coefficients 
corresponding  to  the  greater  local  oriented  energy  are  included  for  reconstruction,  resulting 
in  a  fused  image. 

Figure  6(a)  shows  an  example  of  fusion  of  a  mammographic  image  with  the  modulus  of  the 
dyadic  wavelet  coefficients,  while  Figure  6(b)  presents  an  example  of  fusion  of  a 
mammographic  image  with  the  mammogram  enhanced  by  multiscale  histogram 
equalization.  In  both  cases  the  fused  image  features  enhanced  structures,  while  having  an 
appearance  similar  to  the  original  mammogram. 

2.1.7  Summary 

The  one-dimensional  discrete  dyadic  wavelet  transform  was  extended  to  higher-order 
derivatives  and  even-order  spline  functions,  and  an  improved  initialization  procedure  was 
developed.  In  comparison  to  the  originally  employed  initialization  [7],  our  method  showed 
significantly  better  performance  at  finer  scales  of  analysis  (of  importance  for 
mammographic  features  including  microcalcifications). 
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Figure  6:  (a):  Fused  image  (bottom)  of  an  original  mammogram  (top)  with  the  selected 
modulus  of  dyadic  wavelet  coefficients  (middle),  (b)  Fusion  (bottom)  of  a  dense  mammogram 
(top)  with  a  mammogram  enhanced  by  multiscale  histogram  equalization  (middle). 
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A  multiscale  spline  derivative-based  transform  was  constructed  as  an  approximation  to  a 
steerable  dyadic  wavelet  transform.  The  transform  used  x-y  separable  filters  in  a  perfect 
reconstruction  filter  bank  and  enabled  fast  directional  analysis  of  images  without  the 
introduction  of  artifacts  due  to  translation  and  rotation  invariance.  Such  artifacts  are 
inherent  to  traditional  methods  of  wavelet  analysis. 

Preliminary  image  fusion  results  showed  promise  for  producing  enhanced  images  with  a 
more  familiar  appearance  to  radiologists. 

2.2  A  Parallel  Algorithm  to  Support  Interactive  Wavelet 
Processing  on  a  Radiologist  Workstation 

2.2.1  Introduction 

Diagnosis  of  mammography  is  one  of  medical  imaging  applications  where  digital  image 
processing  is  used  to  extract  either  implicit  or  explicit  characteristics  from  digital  data. 
Many  of  medical  imaging  applications  require  broad  base-level  functionality  for  image 
manipulation,  and  display  to  support  technical  end-users  in  areas  of  Radiology  and  CAD 
(Computer  Aided  Diagnosis).  An  electronic  view-box  for  medical  imaging  analysis  may 
require  advanced  image  processing  capabilities,  including  transformation,  filtering  or 
enhancement  [39].  In  addition,  common  image  manipulation  facilities  such  as  convolution, 
warping,  logical  operations,  comparisons,  and  arithmetic  are  often  desirable  to  exploit  the 
large  dynamic  range  in  contrast  resolution  provided  by  digital  X-ray  detectors. 

In  this  chapter  we  describe  the  design  of  a  high-speed  parallel  algorithm  to  support 
interactive  analysis  via  multiscale  processing.  In  addition  to  writing  multi-threaded  code, 
based  on  kernel  support  for  multiple  threads  of  control  and  changing  the  scheduling  class 
policy  to  real-time  mode,  our  algorithm  required  hardware  and  software  support  for 
high-speed  image  manipulation. 

We  used  a  SPARCstation  20  Model  71  with  an  SX  frame  buffer  to  execute  our 
algorithm.  The  SPARCstation  SX  system’s  direct  accessibility  to  physical  and  display 
memory  provided  flexible  and  powerful  capabilities  for  the  management,  manipulation,  and 
display  of  large  image  matrix  sizes,  typical  of  digital  mammography 

On  top  of  the  SPARCstation  SX  hardware,  foundation  graphics  libraries  were  utilized. 
XGL  and  XIL  accelerate  many  of  graphics  and  X-extension  library  functions.  As  a 
foundation-level  imaging  interface  for  the  SUN  Unix  programming  environment,  the  XIL 
library  defines  how  imaging  operations  such  as  display  and  image  manipulation  are  carried 
out  on  the  workstation  [40]. 

In  this  section,  we  review  the  architecture  of  the  SPARCstation  SX  system,  and 


42 


describe  briefly  the  XIL  foundation  graphics  interface  and  relevant  functions  accelerated  on 
the  SX  system  processor.  We  then  introduce  new  algorithm  to  accomplish  the  forward  and 
inverse  discrete  dyadic  wavelet  transforms  that  exploits  the  parallel  architecture  of  the 
machine. 

2.2.2  Architectural  Overview  of  SPARCstation  SX  System 

The  SPARCstation  SX  graphics  processor  is  a  scalable  graphics  architecture  consisting  of 
an  enhancement  to  memory  controller  chip  on  an  MBus  within  a  SPARCstation.  The 
SPARCstation  SX  system  relies  on  the  notion  that  graphics  operations  become  extensions 
of  native  integer  and  floating-point  operations  of  the  CPU.  The  SPARCstation  SX  system 
employs  a  closely-coupled  architecture,  wherein  a  dedicated  processor  handles  low-level 
display  operations  through  processor-memory  direct  access.  Residing  directly  on  the 
CPU-memory  bus,  the  SPARCstation  SX  memory  controller,  called  the  SMC,  can  perform 
extremely  fast  transfer  of  data  from  physical  memory  (DRAM)  to  the  display  frame  buffer. 
The  close-coupling  of  the  SMC  to  the  SPARCstation  CPU  further  enables  floating 
point-based  computation  to  utilize  the  high-performance  transformation  capabilities  on  the 
SPARC  processor.  Figure  7  illustrates  the  high-level  architecture  of  the  SPARCstation 
MBus  system,  including  CPU,  SMC  processor,  and  video  memory. 


Main  Memory 


The  SMC  processor  lies  on  the  MBus  between  the  SPARC  CPU(s)  and  memory.  In 
addition.  Video  RAM  (VRAM)  is  mapped  into  the  same  address  space,  allowing  the  SMC 
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to  access  and  manipulate  the  display  in  the  same  manner  as  DRAM.  This  addressing  gives 
the  SPARCstation  SX  system  superior  bandwidth  in  moving  large  size  images,  e.g. 
mammograms,  from  memory  to  the  screen. 

Instructions  executing  in  one  or  more  of  the  SPARC  CPUs  make  requests  of  the 
memory  controller  (MC)  portion  of  the  processor.  The  memory  controller  logic  fulfills  CPU 
requests  for  memory  accesses,  such  as  those  required  for  fetching,  loading  and  storing 
instructions. 

When  one  of  the  CPUs  encounters  a  graphics  operation,  it  passes  it  to  the  SMC  pixel 
processor  (SXPP),  which  executes  the  instruction  directly.  When  memory  access  is 
required  to  complete  an  operation,  the  SXPP  makes  a  request  to  the  memory  controller. 

2.2.3  XIL  Imaging  Library 

The  XIL  library  from  Sun  Microsystems  is  a  foundation-level  imaging  interface  providing 
common  functions  required  by  most  imaging  applications.  Such  a  library  is  an  application 
programming  interface  (API),  but  it  is  different  from  “normal”  APIs  because  it  is  hardware 
dependent.  The  XIL  library  contains  three  primary  components:  a  programming  interface 
specification  for  basic  imaging  functionality,  a  high-performance  implementation  of  the 
specification,  and  a  standard  hardware  interface  specification,  which  enables  third-party 
hardware  developers  to  readily  support  the  XIL  library  and  applications  built  on  it  [40]. 
Figure  8  illustrates  the  XIL  library  in  relationship  to  other  foundation- level  interfaces. 


Figure  8:  Solaris  foundation  graphics  libraries  and  layered  interfaces  [40]. 

Monadic  and  dyadic  operations  are  of  particular  interest  within  XIL  library.  Every 
XIL  operator  can  be  considered  an  atom.  Groups  of  XIL  operators  concatenated  in  a 
sequential  manner  are  known  as  molecules.  Molecules  are  important  to  the  SPARCstation 


SX  SMC  processor  because  they  allow  the  device  to  notify  the  XIL  library  at  runtime,  and 
obtain  hardware  support  for  groups  of  chained  operations.  At  runtime,  the  SMC  device 
handler  loads  into  memory  and  notifies  the  XIL  library  of  the  atoms  and  molecules  it 
wishes  to  replace  with  its  own  equivalent  mapping  to  hardware.  The  XIL  library  uses  a 
directed  acyclic  graph  to  hold  a  list  of  these  dependencies  during  execution.  When  the  XIL 
library  encounters  a  known  sequence  of  operations,  it  transfers  the  operation  to  the  SMC 
for  execution.  XIL  applications  can  realize  significant  performance  improvement  from  the 
deferred  execution  [40]. 

2.2.4  An  Algorithm  for  High-Speed  Wavelet  Analysis 

In  frame-based  multi-scale  representations  [32],  spatial  convolution  require  2*  -  1  zeros 
between  each  non-zero  filter  coefficients  in  a  filter  kernel  as  scale  changes.  This  increases 
the  computational  complexity  exponentially.  However,  there  exists  an  algorithm  that 
reduces  the  computational  complexity  by  a  factor  of  2h  The  basic  idea  behind  this 
algorithm  is  to  keep  the  size  of  the  filter  kernels  unchanged  throughout  multiscale  analysis. 
Instead,  we  partition  the  input  image  into  four  sub-regions  so  that  each  region  contains 
only  those  pixels  affected  by  the  filter  kernel  (without  expansion). 

Figure  9  illustrates  differences  between  the  two  approaches:  a  linear  convolution  with 
zero-paddings  in  the  kernel,  and  an  alternative  method  by  splitting  the  convolved  image 
into  four  sub-regions  for  further  convolution  at  finer  levels.  The  former  is  shown  in  Figure 
9(a).  At  level  2,  zeroes  are  inserted  between  the  non-zero  filter  coefficients.  Figure  9(b) 
graphically  illustrates  the  new  method.  First,  the  input  image  is  convolved  with  the 
original  filter  kernel  at  level  1.  The  convolved  image  is  then  divided  into  four  sub-regions 
such  that  each  region  consists  of  only  those  pixels  processed  by  the  convolution  at  the  next 
level  {1-2).  After  this  splitting,  convolutions  are  then  applied  to  each  of  the  sub-regions 
independently.  To  visually  describe  this,  we  show  within  the  input  image  matrix  four 
distinct  pixel  shapes:  dot,  triangle,  square,  and  x,  in  such  a  way  that  pixels  that  belong  to 
the  same  group  are  identified. 

For  example,  let  us  compute  the  value  located  at  coordinate  (2,2)  in  the  input  image 
shown  in  Figure  9(a)  marked  by  a  dot.  All  pixels  within  the  rectangular  area  from  (0,0)  to 
(4,4)  will  participate  in  the  computation.  However,  with  respect  to  traditional  convolution 
the  zeros  in  the  filter  kernel  at  level  2  effectively  mask  out  all  other  pixels  having  shapes 
other  than  a  dot.  Thus  pixels  that  actually  take  part  in  the  computation  are  those  having 
a  dot  shape  located  at  coordinates  (0,0),  (0,2),  (0,4),  (2,0),  (2, 2), (2,4),  (4,0),  (4,2)  and  (4,4). 
The  region  labeled  I  shown  in  Figure  9(b)  at  the  level  2  contains  exactly  those  pixels.  The 
convolution  is  carried  out  separately  in  this  region  and  gives  us  the  same  result  numerically 
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Figure  9;  (a)  Traditional  linear  convolution  (denoted  by  ®)  at  scale  /  =  1,  /  =  2.  One  zero  is 
padded  between  non-zero  filter  coefficients  in  the  filter  kernel  at  scale  2.  (b)  Alternative  approach 
that  avoids  zero-paddings  between  the  non-zero  filter  coefficients.  We  first  convolve  the  input  image 
with  the  original  filter  kernel  at  the  level  1.  The  convolved  image  is  then  divided  into  four  sub- 
regions,  marked  in  the  background  as  I,  II,  III,  and  IV  with  half-grayed  tone,  each  of  which  holds 
pixels  that  are  precisely  affected  by  the  filter  kernel  coefficients  at  the  next  level  I  =  2.  Four  distinct 
pixel  shapes  inside  the  input  image  are  shown:  •,  A,  □  and  x.  Only  pixels  that  have  the  same 
shape  are  processed  by  convolution  at  the  next  level  (denoted  by  ®). 
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as  the  traditional  method. 

Synthesis  is  simply  the  reverse  of  analysis.  In  other  words,  we  first  merge  four 
sub-regions  at  a  coarser  level  into  a  larger  region,  then  convolve  it  with  a  filter  kernel. 
Figure  10  shows  the  synthesis  process. 
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Figure  10:  Multiresolution  block  convolution  for  image  synthesis.  Four  sub-regions  in  the  input 
image  at  a  coarser  level,  marked  as  I,  II,  III,  and  IV,  are  merged  into  a  larger  region  at  a  finer 
level.  Then,  convolution  is  applied  to  the  merged  region. 


For  quantitative  comparison,  let  us  compute  the  computational  complexity  of  each 
method  of  implementation.  For  simplicity,  we  focus  on  the  analysis  only.  Let  the  size  of  the 
input  image  and  the  size  of  the  original  filter  kernel  be  n  and  m,  respectively.  And  let  us 
denote  the  level  of  analysis  by  1.  Then,  in  the  traditional  linear  convolution  case,  the 
computational  complexity  Ca  is 

Ca  —  -l-  2^  —  2)^, 

and  in  our  approach,  the  computational  complexity  Cb  is 

Cb  =  V?  *  2^“^  *  m^. 

Normally,  the  size  of  the  filter  kernel  is  small,  so  we  can  assume  m  is  a  constant.  The 
speedup  by  the  new  approach  is  then 


Algorithm  1  computes  the  two-dimensional  discrete  wavelet  transform  by  this  more 
efficient  method.  The  variables,  n  and  m,  denote  the  size  of  the  input  image  and  the  size  of 
the  filter  kernel,  respectively.  [Filter Kernel]^  y  i®  created  by  the  tensor  product  of  two  ID 
filter  kernels:  X,  Y.  Thus, 

[FilterKernel]^  Y  =  X'  Y. 

defines  a  ROI  of  size  wxh  starting  at  [r  x  tw,  c  x  fi] .  There  are  three  global 
variables  in  the  algorithm:  LEVEL,  WIDTH,  HEIGHT.  By  LEVEL,  we  denote  the  number 
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Algorithm  1  Two-dimensional  discrete  dyadic  wavelet  transform  using  block  splitting 
method _ 

Require:  n  ;  dyadic  number  (n  >  m) 

for  ^  =  0  to  /  <  LEVEL  do 

W  =  WIDTH  of  ROI  } 

h  =  ;  height  of  ROI  } 

for  r  =  0  to  r  <  2^  do 
for  c  =  0  to  c  <  2^  do 

[ROI]^’^  <=  Obtain  a  region  of  interest  from  the  input  image  with  size  of 
w  X  h  starting  at  [r  x  w,c  x  h] 

<=  Split  ([BO/l”;‘  »  [Filter Kernel]^  a) 

Mf]rc  ^  Split  (ifiO/)”;*  .  [FilterKernel]^!) 

[ROI]"*  4=  Split  *  [FilterKernel]„  „') 

{The  process  of  block  splitting  is  shown  in  Figure  9(b).} 

end  for 
end  for 
end  for 


Algorithm  2  Two-dimensional  inverse  discrete  dyadic  wavelet  transform  by  merging  sub- 
regions _ 

Require:  n  :  dyadic  number  (n  >  m) 

for  I  =  LEVEL  —  1  downto  Z  =  0  do 
W  =  WIDTH  ROI  } 

h  =  {h  :  height  of  ROI } 

for  r  =  0  to  r  <  2^  do 
for  c  =  0  to  c  <  2*  do 

[ROI]^’^  ^  Merge  sub-regions  of  [ROI]^’^  *  [Filter Kernel]*fjjj 
[W^if]'^’^  <=  Merge  sub-regions  of  [^2‘ATc  *  [-f'*^^6rAerne/]^,3 
[W^if]'^’^  <=  Merge  sub-regions  of  *  [Filter Kernel]*Qj^ 

[Ror;:*  4=  [ROI]",*  +  +  [wi,f]:f 

{The  process  of  merging  of  sub-regions  is  shown  in  Figure  10.} 

end  for 
end  for 
end  for 
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of  levels  for  the  transform,  and  by  WIDTH  and  HEIGHT  the  width  and  height  of  the 
input  image,  respectively. 

Algorithm  2  describes  our  algorithm  for  computing  two-dimensional  inverse  discrete 
dyadic  wavelet  transform  using  the  method  shown  in  Figure  10.  [Filter Kernel]-)^  Y  is  the 
complex  conjugate  of  [Filter Kernel]^  y- 

We  show  visually  in  Figure  11  the  process  of  a  two  dimensional  forward  wavelet 
transform  using  the  method  described  in  Algorithm  1.  In  Figure  11(b),  we  display  a 
series  of  DC  components  at  levels  /  =  1,  /  =  2  and  I  —  3.  Finally,  Figure  11(c)  shows  an 
enhanced  version  of  the  selected  ROI  shown  in  Figure  11(a)  accomplished  by  execution  on 
the  SPARCstation. 

We  show  in  Table  4  the  computational  costs  in  time  (second)  for  enhancement 
between  methods  of  convolution  by  zero-padding  between  non-zero  filter  coefficients, 
forward  and  inverse  FFT,  and  our  new  multiscale  approach  of  splitting  and  merging 
sub-regions.  Our  sample  input  ROI  was  256  x  256  and  was  computed  up  to  5  levels  of 
analysis.  The  size  of  the  original  kernels  used  for  the  decomposition  was  3  x  3  for  both 
horizontal  and  vertical  directional  filters.  In  the  reconstruction  case,  the  kernels  were  3x5 
in  the  horizontal  direction  and  5  x  3  in  the  vertical  direction. 


level 

Multiscale  Block 

Traditional  convolution 

FFT  approach 

1 

1.18 

0.83 

13.75 

2 

1.81 

2.87 

27.46 

3 

3.71 

7.87 

41.05 

4 

6.29 

26.01 

54.65 

5 

13.03 

93.93 

68.62 

Table  4:  Computational  costs  in  time  (seconds)  of  splitting  (in  decomposition)  and  merging  (in 
reconstruction)  schemes  versus  traditional  linear  convolution  with  zero-padding  between  each  non¬ 
zero  filter  coefficients  and  FFT  implementations. 

The  costs  shown  in  this  table  were  measured  on  a  SPARCstation  20  Model  71  from 
Sun  Microsystems  having  a  75MHz  SuperSPARC-II  processor  with  1MB  SuperCache, 
96MB  RAM,  and  4MB  (60ns)  SX  frame  buffer.  This  table  clearly  shows  that  our  approach 
is  more  efficient  than  existing  convolution  schemes  and  exceeded  an  FFT  by  a  factor  of  5 
for  five  levels  of  processing. 

We  have  developed  a  user  interface,  called  XEnh,  which  provides  an  environment  for 
our  algorithm  using  XIL  foundation-level  imaging  libraries.  This  program  allows  users  to 
select  a  ROI  in  the  main  window  and  carry  out  a  multiscale  contrast  enhancement 
algorithm  (described  in  Chapter  2.1)  within  this  ROI.  Users  are  able  to  set  parameters 
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Figure  11:  (a)  Selected  ROI  containing  an  ill-defined  mass,  (b)  DC  components  at  levels  I  = 
I  =  2  and  I  =  3  after  applying  Algorithm  1.  (c)  Processed  ROI  with  suspicious  area  enhanced 
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This  page  will  provide  you  a  way  to  get  acquainted  with  a  wavelet-based  enliancement 
technique  for  processing  of  manuno grams. 

In  the  Method  section,  our  approach  for  conti  ast  enhancement  of  digital  mammograms  will 
be  explained.  Next,  tlie  Hands-On  section  will  give  you  the  opportunity  to  iry  the 
enliancement  algoritlim  on  preselected  areas  of  mammogiams  m  an  interactive  fashion.  The 
subsequent  Quiz  section  w^iU  enable  you  to  test  the  power  of  enhancement  on  malignant 
mammograms,  and  to  compare  youi'  findings  with  the  actual  diagnosis.  Finally,  tlie 
Conclusions  section  will  summarize  what  we  learned  today. 
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Figure  13:  RSNA  ’96  Exhibition  Home-page  (http://www.iprg.cise.ufl.edu/demo). 
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Level  =  5,  Gain  =  15,  Threshold  =  0.02 


Enhanced 


Original 


interactively,  such  as  the  level  of  transform,  and  values  of  gain  and  threshold.  The  GUI  for 
XEnh  is  shown  in  Figure  12.  The  window  in  the  upper  right  corner  of  the  application  is 
the  area  where  a  selected  ROI  and  its  enhanced  result  are  displayed. 

A  web-site  (http;//www. iprg.cise.ufl.edu/demo)  shown  in  Figure  13  has  been 
established,  and  was  originally  created  for  the  RSNA  (Radiological  Society  of  North 
America)  1996  exhibition  in  Chicago,  Illinois.  In  this  page,  we  provide  a  wavelet-based 
enhancement  technique  for  processing  of  mammograms.  Sections  on  this  page  include 
“Hands-on”  and  “Quiz” ,  where  users  can  process  preselected  or  interactively  selected  ROIs 
of  mammograms  and  obtain  enhanced  results  based  on  the  algorithm  described  in  Chapter 
2.1. 

Figure  14  shows  the  interface  during  the  process  of  enhancement  for  a  sample  case: 
M019  RCC,  where  biopsy  proven  masses  exist.  Figure  14(a)  shows  the  case  with  a 
preselected  ROI.  Figure  14(b)  displays  the  result  of  enhancement  on  the  preselected  ROI 
shown  in  Figure  14(a). 

2.2.5  Summary 

We  have  described  a  parallel  algorithm  to  accomplish  high-speed  interactive  multiscale 
processing  for  enhancement  within  ROIs  in  digital  mammograms.  The  implementation 
relied  upon  foundation-level  libraries,  XIL,  on  top  of  a  SPARCstation’s  SX  frame  buffer. 
We  showed  the  amount  of  speedup  for  the  method  compared  to  traditional  linear 
convolution  and  a  conventional  FFT  approach.  Our  multiscale  approach  employing 
splitting  in  decomposition  and  merging  in  reconstruction  was  shown  to  be  efficient. 

2.3  A  Continuous  Scale  Discrete  Wavelet  Transform  for  the 
Detection  of  Spicular  Masses  in  Mammograms 

2.3.1  Introduction 

The  Continuous  Wavelet  Transform  (CWT)  has  been  shown  to  be  an  effective  tool  for  the 
analysis  of  non-stationary  signals  [8]  [41].  A  CWT  performs  a  time-scale  analysis,  zooms  in 
on  singularities,  captures  slow  variations  and  provides  a  scale-invariant  interpretation  of 
shapes. 

The  importance  of  scale  has  been  acknowledged  in  computer  vision  for  a  long  time. 
For  example,  in  1980,  Marr  and  Hildreth  developed  a  multiscale  edge  detection  algorithm 
modeled  on  the  human  visual  system  [42].  More  recently,  it  has  been  shown  that  wavelet 
theory  can  provide  a  solid  foundation  for  multiscale  analysis.  For  example,  tracking 
wavelet  transform  maxima  across  scales  can  allow  us  to  discriminate  intensity  profiles  of 
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distinct  types  of  edges. 

In  image  processing,  Discrete  Wavelet  Transforms  (DWT)  are  more  popular  than 
CWTs  because  most  signals  are  digitally  sampled.  Mallat  first  developed  a  fast  DWT 
algorithm  based  on  multi-resolution  analysis  [41].  It  computed  the  wavelet  coefficients  on  a 
sampled  grid  (a, b)  defined  hj  a  2^k,  where  i,keZ.  Its  wavelet  representation  is 

complete  and  concise.  The  *^a  trous”  algorithm  calculates  wavelet  coefficients  at  a  denser 
sample  grid,  a  =  2\b  =  k  [13].  An  important  advantage  of  the  “a  trous”  algorithm  for 
image  analysis  is  that  it  is  shift-invariant. 

Visual  features  occur  at  distinct  scales.  For  feature  analysis,  it  is  desirable  to  detect 
features  at  an  arbitrary  scale.  M.  Unser  et.  al.  [43]  [44]  have  proposed  several  fast  and 
efficient  algorithms  for  calculating  wavelet  coefficients  at  arbitrary  scales.  In  these 
algorithms,  Q  scales  were  inserted  into  each  octave  to  obtain  a  denser  sampling  of  scale. 
The  wavelets  at  each  of  the  Q  scales  in  the  first  octave  were  approximated  by  their 
projections  in  scale  space,  and  called  auxiliary  wavelets.  Since  the  auxiliary  wavelets  lie  in 
scale  space,  auxiliary  wavelets  at  the  next  higher  octave  were  obtained  by  using  the 
standard  two-scale  relation. 

In  this  chapter,  we  introduce  a  special  type  of  mother  wavelet  that  allows  for  an 
arbitrary  sampling  of  the  scale  parameter.  We  discuss  how  to  create  such  a  wavelet,  which 
is  an  approximation  of  a  known  mother  wavelet.  We  also  develop  a  method  to  reconstruct 
an  original  image  from  its  wavelet  coefficients  generated  from  an  arbitrary  scale  analysis. 
We  demonstrate  the  importance  of  arbitrary  scale  analysis  by  applying  the  method  to 
problems  of  detecting  breast  mass  lesions  of  arbitrary  size  in  digitized  mammograms. 

The  chapter  is  organized  as  follows.  In  Section  2.3.3,  we  first  derive  the  Continuous 
Scale  Discrete  Wavelet  Decomposition  from  the  CWT.  In  Section  2.3.4,  we  show  how  to 
reconstruct  a  signal  from  its  wavelet  coefficients  at  an  arbitrary  scale  a  and  at  scales  2“*a, 
where  i  G  Z,  2“*o  >  Oq  and  oq  is  some  initial  scale  of  analysis.  In  Section  2.3.5,  we  discuss 
an  algorithm  using  wavelets  that  are  dth  derivatives  of  B-spline  wavelets  [7]  [31].  We  then 
define  a  general  form  of  filters  to  implement  our  multiscale  analysis  and  synthesis  in  the 
frequency  domain.  In  Section  2.3.6,  we  show  the  results  of  applying  our  algorithm  to 
difficult  mammography  cases  of  breast  cancer.  The  ability  to  find  and  process  features  at 
an  arbitrary  scale  is  shown  to  have  significant  advantages  over  standard  dyadic  algorithms. 
In  one  case,  the  algorithm  detected  the  center  of  a  mass  that  could  not  be  seen  using 
conventional  window  and  leveling  techniques  with  12-bit  contrast  resolution.  In  Section 
2.3.7,  we  discuss  the  relation  of  our  decomposition  scheme  with  the  oblique  projection 
scheme  proposed  by  M.  Vrhel,  M.  Unser  et  al.  [43]  and  show  that  they  are  identical  for  a 
limiting  case.  Finally,  in  Section  2.3.8,  we  summarize  our  results. 
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2.3.2  Definitions  and  Notation 

In  this  chapter,  a  letter  in  lower  case  represents  a  function  or  a  sequence  in  the  time  or 
spatial  domain.  Its  upper  case  represents  the  corresponding  frequency  response,  i.e. 
p{k)  -H-  P{oj),  •0(t)  -H-  ^(w).  ■0a(t)  =  '0(|)  •H-  Given  a  sequence  x{n),  a 

function  is  reconstructed  from  x{n)  using  an  ideal  lowpass  filter. 

The  symbol  will  be  used  for  three  distinct  types  of  convolutions:  linear 
convolution,  mixed  convolution  and  discrete  convolution.  For  functions  /  and  g  defined  on 
TZ,  denotes  traditional  convolution 

/OO 

fiT)g{t  -  T)dT. 

■OO 

The  mixed  convolution  between  a  sequence  b{k)f.^2:  ^  function  /  defined  on  71  is 

OO 

((-•/)(*)=  E  mnt-k). 

k=—oo 

The  discrete  convolution  between  two  sequences  a  and  b  is  defined  as 

OO 

{a*b){n)—  ^  a{k)b{n  —  k). 

k—  —  oo 

Whenever  it  exits,  the  convolution  inverse  b~^  of  a  sequence  b  is  defined  to  be 

*b)  =  5{k). 


Function  is  Schoenberg’s  central  B-spline  of  degree  n,  which  is  obtained  from  the 
(n  +  l)-fold  convolution  of  a  unit  rectangular  pulse: 


t  + 


n  +  1 
2 


t  G  7^, 


where  [t]+  ==  max{0,t}.  The  Fourier  transform  of  is  given  by  /7"'(a;)  == 

Discrete  B-spline  sequences  are  obtained  by  sampling  the  continuous  B-spline  functions, 
iov  k  e  Z  .  Note  is  stable  for  any  value  of  n  [45]. 

Finally,  the  function  Rect{u))  is  defined  as 


Rect{ix) 


f  1 

when  |a;| 

1  0 

when  jo; 
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2.3.3  A  Continuous  Scale  Discrete  Wavelet  Transform 

A  wavelet  transform  is  a  linear  operation  that  projects  a  signal  onto  a  set  of  basis 
functions,  which  are  dilated  versions  of  some  mother  wavelet.  The  mother  wavelet  is  a 
function  ^  e  that  satisfies  <  +oo.  This  requires  that  the  wavelet  has 

sufficient  decay,  and  that  =  0.  The  continuous  wavelet  transform  of  a  function 

at  scale  a  and  shifting  parameter  b  is  defined  as 

Wf{a,b)  =  l^ 

Of  course,  it  has  been  shown  that  a  wavelet  transform  satisfies  energy  conservation  and  can 
be  perfectly  reconstructed  from  its  wavelet  transform  coefficients  [8]  [46].  The  Continuous 
Scale  Discrete  Wavelet  Transform  (CSWDT)  is  defined  as  the  CWT  sampled  along  the 
shifting  parameter  b, 

Wfia,n)^J  (60) 

In  practice,  data  is  acquired  digitally.  Furthermore,  we  assume  that  a  signal  / (t)  is 
sampled  at  its  Nyquist  rate  fs  so  that  it  may  be  recovered  from  its  discrete  samples.  Since 
f{t)  is  band-limited,  W/(a,  b)  is  also  band-limited,  and  can  be  recovered  from  its  samples 
W/(a,n). 


f(t)  - ► 


Sk5(t-k) 


Figure  15:  Data  acquisition  processing. 

Processing  for  data  acquisition  is  shown  in  Figure  15  (^(t)  is  an  anti-aliasing  filter). 
Since  the  sampled  signal  x(n)  is  of  finite  resolution,  this  puts  a  limit  on  the  finest  scale  that 
can  be  computed.  Let  us  introduce  a  compactly  supported  scaling  function  that 
satisfies  the  following  conditions  [43]  : 

1.  Riesz  basis  condition:  A  <  +  2fc7r)p  <  B,  where  A  and  B  are  strictly 

positive  constants; 

2.  Order  property:  $(0)  =  1,  $^”^^(2A;7r)  =  0,  A:  G  .2,  A:  7^  0,  for  m  =  0, ...,  A"  —  1; 

3.  Two-scale  relation:  (j)  (|)  =  Ylkez  ~  ^)’ 
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4.  Frequency  property:  4>(a;)  ^  0,  for  |a;|  <  tt,  and  $(a;)  should  be  very  small  for 

|a;|  >  TT. 

Constraint  1  implies  that  (f){t  —  k)  is  a  Riesz  basis  of  the  subspace 

v<l>  =  {/W  =  c{k)(l){t  -k)  c  e  I2}  for  f{t)  G 

kez 

and  that  is  a  well  defined  (closed)  subspace  of  C2  [47].  Property  2  implies  that  (t){t) 
reproduces  all  polynomials  of  degree  N-1  [48].  Constraint  3  is  the  well-known  two-scale 
relation.  Condition  4  shows  that  can  be  an  anti-aliasing  filter. 


f(t) 


Ik5(t-k) 


Figure  16:  Initial  conditions. 


Typically,  the  anti-aliasing  filter  in  the  acquisition  device  is  unknown.  For  simplicity, 
we  assume  that 


x{n)  = 


f{t)(j){n  -  t)dt, 


(61) 


which  is  shown  in  Figure  16.  This  is  the  initial  condition  typically  found  in  the  traditional 
DWT  analysis  literature.  When  the  impulse  response  of  the  anti-aliasing  filter  (i.e.  ideal 
low-pass  filter)  is  known,  a  pre-filter  and  a  post-filter  are  used  as  shown  in  Figure  17 
(Please  See  Appendix  A.l  for  details). 


x(n)  - 


Z  ^  *  Wavelet  decomposition 

- - ► 

Z 

I  and  reconstruction 

x(n) 


Figure  17:  Wavelet  processing  with  alternative  initial  conditions. 

The  smoothing  operator  Sa  is  defined  as  Saf{n)  =  jToof{t)<P  (^)  dt.  Then 
Sif{n)  =  x{n),  where  a —  I  is  the  finest  scale  at  analysis.  As  a  increases,  Saf{n)  becomes 
a  coarser  representation  of  f{t).  When  a  goes  to  infinity,  Saf{n)  approaches  a  constant 
equal  to  the  average  value  of  f{t).  In  practice,  the  length  of  x{n)  defines  the  coarsest  scale. 
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Lemma  2.1  Given  x{n),  there  is  one  and  only  one  f{t)  such  that 

x‘^{t)  =  f  *  <p{t).  (62) 

Proof:  Since  X^{uj)  is  band-limited  to  [— 7r,7r]  and  $(a;)  ^  0  for  w  G  [— 7r,7r], 

p(,  ■ 

Theorem  1  Given  x{n)  =  x‘^{t)\t=n>  if 

1.  is  band-limited  between  [— 7r,7r]  when  a  >  Sq; 

2.  a>  So- 

Then,  the  CSDWT  can  be  calculated  by 

Waf{a,  n)  =  J2p{a,  k)x{n  -  k),  (63) 

kez 

where 

P^{a,uj)  =  ^^  (64) 

$(u;) 

and 

/OO 

p^{a,  t)e~^'^^‘*'dt. 

■OO 

Proof:  Note  that  fpa{t)  and  p'^{a,  t)  are  band-limited  at  [-tt,  tt].  And  p{a,  n)  =  p'^{a,  t)\t=:n- 
From  Equation  (64), 

^a(a;)  =  P‘^(a,a;)$(a;)  (65) 


'fpait) 


/: 


u)4>{t  —  u)du. 


(66) 


If  we  substitute  Equation  (66)  into  the  previous  CSDWT  definition  (Equation  (60)),  we 
obtain 


Waf{a,n) 


/  f{t)'ipa{n-t)dt 
J  — OO 

/OO  poo 

f(t)  /  p'^(u,  u)(j){n  —  t  —  u)dudt 

-OO  J  — OO 

poo  poo 

/  (  /  /(^)</’(^  -t  -  u)dt)p‘^{a,  u)d'i 

J —OO  J —OO 

poo 

/  x^{n  —  u)p^{a,  u)du. 


(67) 
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(68) 


Since  x‘^{t)  and  p‘^{a,t)  and  band-limited  to  [— 7r,7r],  we  need  to  show  that 

/OO 

x‘^{n  —  u)p‘^{a,  u)du  =  '^^p{a,  k)x{n  —  k). 

■°°  kez 

If  we  look  in  the  frequency  domain,  we  see  that 

/OO 

x^{n  —  u)p'^{a^u)du}  =  (X^(a;)P'^(a, w))  *  27rh27r(w)- 

-OO 

Since  x'^{t)  andp‘^(a,  t)  are  band-limited  to  [— tt,  tt], 

/OO 

x^{n  —  u)p‘^{a,  u)du}  =  {X‘^{u)  *  27r(52,r(w))(P^(a,  <^)  *  27r(527r(nj)) 

■OO 

=  X{u))P{a,co). 

Therefore,  we  can  claim  that  Equation  (68)  is  indeed  true.  Substituting  Equation  (68)  into 
Equation  (67),  we  obtain 

Waf{a,  n)  =  '^pia,  k)x{n  -  k). 


kez 


Thus,  p(a,  k)  of  length  N  can  be  calculated  as  follows: 

1.  For  uj  G  [— TT,  tt).  Evaluate  P{a,Lo)  =  P‘^{a,Ljj)  =  at  a;  =  ,  where 

—N  <  n  <  N] 


2.  Note  that  P{a,uj)  is  a  27r  periodic  function,  calculate  P{a,uj)  for  co  G  [0,27r); 

3.  Calculate  p(o,  k)  from  the  IFFT  of  P(a,a;). 


Usually,  a  wavelet  'ip{t)  is  not  compact  in  the  frequency  domain.  However,  ^(cu)  is 
local  and  we  can  select  an  sq  such  that  ^so(^)  =  acceptable 

approximation  error  onto  the  scaling  space.  The  approximation  error  is  defined  as  the  ratio 
of  the  energy  of  ^(ow)  outside  [-7r,7r]  to  the  total  energy.  It  decreases  as  Sq  increases. 

Next  we  show  that  our  'ipsoit)  is  also  a  mother  wavelet.  Since  ^p{t)  satisfies  the 
“admissibility  condition” , 


-^\^{uj)\‘^dLO 


/: 


Isowl  ^|\I/(so(^)p(i(sow) 


cx: 

TT 


\u 

\(jj 


^\^{SQC0)\‘^duJ 

~^\^{SQLj)\^duJ 


^\^{So(jj)\‘^dLO 


< 

< 

< 

< 

< 


OO 

OO 

OO 

OO 


OO. 


(69) 
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Th6r6for6,  'ipso  (^)  th6  ^^adinissibility  condition  .  W3.vcl6t  coefficients  nt  arbitrary 

scale  a  =  ssq  can  be  calculated  directly  from  Equation  (63).  However,  wavelet  coefficients 
at  scale  a  =  2*So  can  be  calculated  efficiently  using  the  fast  algorithm  discussed  in  the  next 
section. 


2.3.4  An  Efficient  Algorithm  for  Analysis  and  Synthesis 

When  there  is  knowledge  about  the  best  scale  in  which  to  detect  a  certain  feature,  it  is 
desirable  to  identify  the  feature  at  that  scale,  enhance  the  feature  and  reconstruct  the 
signal.  We  now  derive  a  computationally  efficient  scheme  to  decompose  a  signal  onto  scale 
(  a  =  2^ So  )  and  reconstruct  the  signal  at  that  scale. 

Theorem  2  Let  ipaoit)  =  ~ 

L-l 

fc=0 

Proof:  Since 

$(soa;)  =  F(so,a;)$(a;),  (70) 


we  can  substitute  u)  with  2^lo, 

^{2^ Sou;)  =  P{so,  2^u)^2^uj).  (71) 

Then,  using  the  two  scale  relation  to  simplify  P{so,  2^a;)<l>(2^a;), 

P{so,2^to)^2^uj)  =  P{so,2^uj)H{2^-^u;)^2^-^u) 

L-l 

=  P(so,2^u;)(Y[H{2^u;))^u;).  (72) 

k=0 

Substituting  Equation  (72)  into  Equation  (71), 

L-l 

^{2^ Sou;)  =  P^so,  2^u;){ll  H{2’^u;))^u;).  (73) 

A:=0 


If  we  let  K{u)  be  a  27r  periodic  function.  For  u;  e  [-tt,  tt] 

l-\Hiu;)f 


K{u)^ 


P(so,uj) 


(74) 
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Figure  18:  Filter  bank  for  CSDWT  decomposition  and  reconstruction  for  3  levels  of  analysis. 


Thus  a  signal  can  be  decomposed  and  reconstructed  at  an  arbitrary  scale  as  shown  in 


Figure  18.  Note  that  the  approximation  error  of  is  dependent  upon 


J+°° 

5'(5ow)da; 


and 


.  The  amount  of  error  can  be  reduced  by  judicious  selection  of  sq  and  $(t). 


For  the  dyadic  wavelet  transform,  we  have  T(2a;)  =  G(a;)$(a;).  Therefore 


G{io) 


In  the  CSDWT  case,  we  simply  write 

P{so,u)Rect  = 


T(2a;) 

«I>(a;)  ■ 

(75) 

Vtt/ 

(76) 

Thus  when  sq  =  2,  G{u)  =  P{so,u)  for  tt  G  [-7r,7r].  Since  G{lo)  is  27r  periodic, 

G{u)  =  P{so,oj)  for  all  oj.  In  the  next  section,  we  provide  some  examples  for  the  design  of 
CSDWT  filters. 


2.3.5  Examples  of  CSDWT  Filters 


Spline  functions  have  played  a  significant  role  in  wavelet  theory.  The  well-known  Haar 
transform  corresponds  to  a  spline  of  order  n  =  0.  Battle  and  Lemarie  independently 
constructed  orthogonal  spline  wavelets  using  symmetric  basis  functions  [49].  Chui  and 
Wang  introduced  the  B-spline  wavelets  of  compact  support  [50].  Unser  et.  al.  first 
developed  shift-orthogonal  wavelets  using  splines  [51].  In  this  section,  we  derive  the  filters 
for  a  continuous  scale  DWT  using  the  wavelets  originally  suggested  by  Mallat  in  1992  [7] 
and  extended  by  Laine  and  Koren  in  1996  [31]. 

The  scaling  function  for  a  continuous  scale  transform  is  defined  as 


$(a;)  = 


(77) 
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Its  wavelet  is  written  as 


(78) 

and 

(79) 

The  two-scale  relation  becomes 

(80) 

and 

P(a;)  =  e^^°‘^cos^  , 

(81) 

where  to  =  |  when  n  is  odd  and  level  L  =  0.  Otherwise,  to  =  0.  t^(t)  is  a  d-th  derivative  of 
a  spline  of  order  n  d  —  1  and  4>{t)  is  a  spline  function  of  order  n  1.  The  e  ^  °  term  is 
introduced  because  central  B-splines  have  knots  at  i  e  for  n  even  and  at  z  -1- 1  for  n  odd. 
In  our  fast  decomposition  and  reconstruction  algorithm,  for  oj  G  [-7r,7r],  we  can  write 

Ad2n+2d  sin^+‘^  (^) 

(82) 

H{co)  =  (1) 

(83) 

K{uj)  =  jj,  X  • 

P{sq,u) 

(84) 

Where,  h  =  ^  when  d  is  odd  and  level  L  =  0,  and  ti=0  otherwise.  Notice  that  when 
So  =  2,  it  is  easy  to  show  that  the  filters  P(so,u;),  H{(^)  and  K{u))  are  the  same  as  the 
traditional  DWT  filters  G{uj),  H{uj)  and  K{lo). 

As  an  example,  we  constructed  and  applied  a  spline  wavelet  with  d  =  2  and  n  =  3  to 
decompose  a  sample  of  white  noise  onto  scale  14.25.  We  then  reconstructed  the  signal  from 
the  same  scale.  Multiscale  coefficients  of  the  decomposition  are  shown  in  Figure  19. 
Remarkably,  the  reconstruction  error  was  less  than  10“^^,  which  was  due  principally  to  the 
floating  point  accuracy  of  our  digital  computer,  a  Sun  SparcStation  20. 

In  the  next  section,  we  apply  this  algorithm  to  a  digital  radiograph  of  a  phantom  and 
sample  digitized  mammogram  cases. 
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Fast  decomposition  onto  scale  14.25 


50  100  150  200  250  300  350  400  450  500 


Figure  19:  Sample  of  fast  decomposition  algorithm  of  white  noise  at  scale  14.25. 


2.3.6  Application  of  the  CSDWT  Algorithm 

X-ray  images  of  an  RMI  model  156  phantom  (Radiation  Measurements  Inc.,  Middleton, 
WI)  were  investigated.  The  RMI  phantom  contains  five  masses  and  six  fibrils  which  mimic 
two  objects  of  interest  in  mammography.  The  radiographic  image  quality  of  this  phantom, 
as  well  as  the  total  number  of  objects  which  are  deemed  to  be  visible,  form  a  key 
component  of  the  mammography  accreditation  program  administered  by  the  American 
College  of  Radiology  (ACR).  The  RMI  phantom  is  designed  so  that  at  least  one  object  (i.e. 
smallest  mass)  in  each  category  is  not  visible,  and  one  object  is  at  the  borderline  of  the 
visibility  threshold  [52] . 

Radiographic  images  of  the  RMI  phantom  were  obtained  using  a  Digital  Spot 
Mammography  (DSM)  system  attached  to  a  Lorad  Breast  Biopsy  System  (Lorad 
Corporation,  Danbury,  CT).  Due  to  the  limit  of  the  detector  area,  only  regions  containing 
masses  at  the  right-lower  corner  were  captured  (with  technique  factors  of  22  kVp  and  112 
mAs)  as  shown  in  Figure  20(a).  The  smallest  mass  in  the  schematic  representation  of  the 
insert  of  the  phantom  (Figure  20(e))  can  not  be  seen  using  conventional  window  and 
leveling  at  12-bit  contrast  resolution. 

In  this  study,  the  wavelet  applied  was  a  second  derivative  of  a  spline  of  order  5  (n=4. 
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d=2).  The  scaling  function  was  a  cubic  spline.  The  2-D  wavelet  filters  were  derived  from 
the  1-D  design  using  a  rotated  circularly  symmetric  window  [53].  The  minimums  of  wavelet 
coefficients  computed  at  scale  53.8  are  shown  in  Figure  20(b).  A  soft  threshold  0.01  was 
then  applied.  Note  that  all  the  masses  were  correctly  detected.  Figures  20(c)  and  20(d) 
show  that  a  dyadic  wavelet  fails  to  detect  all  the  masses  within  the  phantom. 

Mammograms  obtained  from  our  local  database  containing  lesions  of  varying  density 
are  described  in  Table  5. 


Case 

Shape 

Margins 

Density 

Size 

Biopsy 

Difficulty 

rml041 

Irregular 

Indistinct 

Equal 

12mm 

Malignancy 

subtle 

lcc046 

Irregular 

Indistinct 

High 

10mm 

Malignancy 

Mod.  Difficult 

lcc002 

Irregular 

Spiculated 

High 

10mm 

Malignancy 

Difficult 

lml015 

Round 

Obscured 

High 

15mm 

Malignancy 

Difficult 

Table  5:  Description  of  mammograms 


Original  images  rml041  and  lcc046  are  shown  in  Figure  21.  In  each  case,  the  scale 
space  was  systematically  searched  (interactively)  to  find  an  optimum  scale  for  detection. 
Then,  chains  of  wavelet  coefficient  minima  at  the  selected  scale  were  labeled  and  classified 
as  shown  in  Figure  22.  For  proper  labelling,  the  set  of  maxima  were  first  thinned  to  be 
exactly  1  pixel  wide.  The  labelling  process  first  classified  each  maxima  as  one  of  the 
following; 

1.  JUNCTION:  If  a  minimum  had  three  or  more  (8-connected)  neighbors; 

2.  TERMINAL:  If  a  minimum  had  zero  or  one  neighbor,  or  exactly  two  neighbors  at 
least  one  of  which  is  a  JUNCTION; 

3.  EDGE;  If  a  minimum  had  exactly  two  neighbors,  neither  a  JUNCTION. 

A  set  of  contour-level  labels  were  generated,  using  the  TERMINALS  to  find  endpoints. 
When  all  TERMINALS  had  been  processed,  any  EDGEs  not  labelled  were  linked  into 
contour-level  labels.  The  JUNCTIONS  were  then  processed  to  produce  group-level  labels, 
each  representing  a  set  of  8-connected  chains.  Each  contour  was  examined  to  determine  if 
it  should  be  split  to  form  chain-level  labels.  Finally,  local  standard  deviation  of  the  gray 
level  intensity  values  for  all  pixels  in  each  contour  was  then  calculated.  Each  chain  was 
assigned  to  a  Group  as  shown  in  Figure  23.  Note  that  chains  projecting  from  the  center  of 
the  mass,  follow  the  spicular  distortion  near  the  mass.  In  each  case,  the  algorithm  detected 
each  mass,  at  scale  124.8  and  scale  112,  respectively.  However,  for  rml041,  the  detection 
algorithm  also  picked  up  an  object  in  the  right-lower  Conner,  which  looked  quite  similar  to 
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Figure  20;  (a)  X-ray  image  of  the  RMI156  phantom;  (b)  Coefficient  minima  obtained  at 
scale  53.8;  (c)  Coefficient  minima  obtained  at  scale  32;  (d)  Coefficient  minima  obtained  at 
scale  64;  (e)  Schematic  representation  of  mammographic  features  within  the  phantom. 
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Figure  22:  (a)  Labeling  of  wavelet  coefficient  minima  within  selected  ROI  from  rml041  (soft- 
thresholding  applied).  JUNCTIONS"  are  shown  in  BLUE,  ‘^TERMINALS"  in  GREEN, 
and  “CONTOURS”  in  red;  (b)  Labeling  of  wavelet  coefficient  minima  within  ROI  from 
lcc046;  (c)  4X  magnification  of  the  center  cluster  representing  a  mass  shown  in  (a);  (d)  4X 
magnification  of  the  center  cluster  representing  a  mass  shown  in  (b). 
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(c)  (d) 

Figure  23:  (a)  Wavelet  coefficients  of  rml041  at  scale  124.8;  (b)  Color  coding  of  a  local 
texture  measure  computed  normal  to  the  contour  of  a  spicular  distortion  around  the  mass 
detected  in  rml041.  Note  consistency  of  measure  for  contours  identified  within  each  cluster; 
(c)  Wavelet  coefficients  of  lcc046  at  scale  112;  (d)  Color  coding  of  a  local  texture  measure 
computed  normal  to  the  contour  of  a  spicular  distortion  around  the  mass  shown  in  lcc046. 


69 


Figure  25:  (a)  Labeling  of  wavelet  coefficient  minima  within  selected  ROI  from  lcc002  (soft- 
thresholding  applied).  junctions”  are  shown  in  BLUE,  ‘‘‘'TERMINALS”  in  GREEN,  and 
“CONTOURS”  in  red;  (b)  Labeling  of  wavelet  coefficient  minima  within  ROI  from  lml015; 
(c)  4X  magnification  of  the  center  cluster  representing  a  subtle  mass  shown  in  (a);  (d)  4X 
magnification  of  the  center  cluster  representing  a  subtle  mass  shown  in  (b). 


(c)  (d) 

Figure  26:  (a)  Wavelet  coefficients  of  lcc002  at  scale  126.4;  (b)  Color  coding  of  a  local  texture 
measure  computed  normal  to  the  contour  of  a  spicular  distortion  around  a  mass  shown  in 
lcc002.  Note  consistency  of  measure  for  contours  identified  within  each  cluster;  (c)  Wavelet 
coefficients  of  lml015  at  scale  89.6;  (d)  Color  coding  of  a  local  texture  measure  computed 
normal  to  the  contour  of  a  spicular  distortion  around  the  mass  shown  in  lml015. 


a  mass.  Additional  texture  information  was  needed  to  separated  it  from  the  true  mass. 
Cases  lcc002  and  lml015  are  shown  in  Figure  24.  These  mammograms  were  rated  as 
“difficult”  cases  by  an  expert  mammographer.  A  selected  ROI  within  each  mammogram 
was  decomposed  onto  an  optimum  scale  126.4  and  89.6  respectively.  Their  wavelet  minima 
are  shown  graphically  in  Figure  25.  Detection  results  are  shown  in  Figure  26.  Both  masses 
were  detected.  It  is  important  to  note  that  in  all  cases,  the  optimum  scale  found  to  detect 
the  mass  was  not  power  of  2.  For  comparison  with  traditional  methods,  Figure  27  shows 
failure  to  detect  the  mass  at  dyadic  scales.  In  addition,  tests  using  the  RMI  phantom 


(a)  (b) 

Figure  27;  (a)  Wavelet  coefficients  for  lml015  at  scale  64;  (b)  Wavelet  minima  in  lml015. 

showed  that  the  detection  algorithm  is  quite  sensitive.  Additional  results  on  real 
mammograms  appear  very  promising.  Masses  within  dense  mammograms  of  subtle 
visibility  were  detected! 

2.3.7  Frequency  Approximation  vs  Projection  Methods 

In  this  section,  we  discuss  the  relation  between  this  algorithm  with  the  oblique  projection 
method  proposed  by  Unser  et.  al.  etc.  [43].  In  addition,  we  show  in  Theorem  3  that  they 
are  equivalent  in  an  extreme  case.  For  a  scaling  function  (j),  we  choose  an  analysis  function 
^  (pi.  Next  the  dilated  wavelet  is  projected  into  y((/>)  such  that  the  approximation 
error  is  orthogonal  to  V{4>i).  Here 

Pa{k)  =  {qa*Qi2)ik),  (85) 

^An  analysis  function  is  defined  as  a  function  which  satisfies  only  scaling  function  properties  (1)  and  (2), 
as  described  in  section2.3.3  at  this  chapter. 
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Qaik)  =<  'ifait),  -  k)  >, 

(86) 

ai{k)  =<  f)i{t  -  k),(t){t)  >, 

(87) 

qi2  = 

(88) 

When  G  Pa  is  an  orthogonal  projection.  Ideally,  (j)i{t)  should  be  chosen  such 

that  the  oblique  projection’s  approximation  error  is  very  close  to  that  of  an  orthogonal 
projection,  and  the  computation  is  simple  and  fast.  We  chose  tpaii)  =  i’  is 

selected  based  on  some  estimate  of  an  acceptable  projection  error.  The  parameter  5o  sets 
the  finest  scale  resolution,  and  the  continuous  scale  discrete  wavelet  transform  of  f{t)  is 

yVaf{a,n)^'^Pa{k)x{n-k).  (89) 

kez 

Theorem  3  Let  4>i{t)  =  When  n-\-m  +oo,  the  oblique 

projection  method  is  equivalent  to  the  frequency  approximation  method. 

Please  see  proof  in  Appendix  A. 2.  As  claimed  in  Theorem  3,  the  frequency 
approximation  method  is  equivalent  to  the  projection  method  in  one  limiting  case.  The 
advantage  of  the  oblique  projection  method  is  that  it  can  be  implemented  by  FIR  and  HR 
filters  to  achieve  0(N)  complexity  per  scale.  However,  the  frequency  based  approximation 
method  is  easy  to  expand  for  multi-dimension  signal  analysis.  Choosing  a  higher  order 
scaling  function  does  not  afi'ect  the  computation  complexity  adversely. 

2.3.8  Summary 

We  showed  that  the  continuous  wavelet  transform  of  a  band-limited  signal  can  be 
estimated  from  samples  of  a  signal.  If  a  mother  wavelet  is  compact  in  the  frequency 
domain,  effectively  no  approximation  error  was  observed  experimentally. 

We  described  a  modified  dyadic  wavelet  transform,  which  was  able  to  decompose  a 
signal  onto  an  arbitrary  scale  and  reconstruct  perfectly  the  signal  from  that  scale.  We 
presented  a  special  class  of  spline  wavelet  filters,  which  allowed  us  to  use  a  high  order 
scaling  function  without  significantly  increasing  the  computational  complexity  of  analysis. 

The  importance  of  arbitrary  scale  analysis  was  demonstrated  by  digital  radiographs  of 
a  mammography  phantom  and  digitized  mammograms.  The  study  showed  that  the 
algorithm  was  able  to  detect  very  subtle  masses,  which  were  rated  to  be  almost  invisible  by 
radiologist  specializing  in  mammography. 
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2.4  Enhancement  of  Mammograms  from  Oriented  Information 

2.4.1  Introduction 

Mammograms  can  detect  tumors  that  are  an  eighth  of  an  inch  in  diameter,  while  manual 
examination  usually  fails  to  detect  tumors  smaller  than  a  half-inch.  Reliable  diagnosis  by 
radiographs  of  malignant  breast  disease  depends  on  observing  local  and  distant  changes  in 
tissues  produced  by  disease.  Of  the  visual  signs  of  cancer  found  by  radiologists,  the 
primary  signs  of  masses  and  calcifications  are  most  important  [54,  55].  Unfortunately,  at 
the  early  stages  of  breast  cancer,  these  signs  are  very  subtle  and  varied  in  appearance, 
making  diagnosis  difficult  and  challenging  even  to  specialists  [54,  56]. 

The  primary  radiographic  signs  of  cancer  are  related  to  tumor  mass  and  its  density, 
size,  shape,  borders  and  calcification  content.  Extraction  of  these  features  and 
enhancement  of  them  may  assist  radiologists  to  locate  suspicious  areas  more  reliably  [57]. 

Multiscale  representations  based  on  wavelets  have  been  carried  out  for  mammographic 
feature  analysis  [32,  58,  59].  Laine  et  al.  [32]  used  two  overcomplete  multiscale 
representations  for  contrast  enhancement.  Mammograms  were  reconstructed  from 
transform  coefficients  modified  at  each  level  by  nonlinear  weighting  functions.  Qian  et  al. 
[58]  introduced  tree-structured  nonlinear  filters  for  micro  calcification  cluster  detection.  An 
image  was  enhanced  by  tree-structured  nonlinear  filters  with  fixed  parameters  and  adaptive 
order  statistic  filters.  Richardson  Jr.  [59]  applied  linear  and  nonlinear  filtering  approaches 
to  the  analysis  of  mammograms.  Here,  a  linear  multiscale  decomposition  was  obtained  via 
a  wavelet  transform;  a  nonlinear  multiscale  decomposition  employed  a  “mean  curvature 
partial  differential  equation”  filter  and  “weighted  majority-minimum  range”  filter.  In 
addition,  Li  et  al.  [60]  extended  a  conventional  multiresolution  wavelet  transform  into  a 
multiresolution  and  multiorientation  wavelet  transform.  They  applied  directional  wavelet 
analysis  to  capture  orientation  information  within  each  mammogram. 

Freeman  and  Adelson  [20]  first  proposed  the  concept  of  steerable  filters  and  applied  it 
to  several  problems  in  the  area  of  computer  vision.  With  a  set  of  “basis  filters”,  one  can 
adaptively  steer  a  filter  along  any  orientation.  Hilbert  transform  pairs  were  constructed  to 
find  a  local  “oriented  energy”  measure  and  dominant  orientation. 

Kass  and  Witkin  [61]  developed  an  algorithm  for  estimating  the  orientation  of  texture 
patterns.  An  orientation  pattern  was  decomposed  into  a  fiow  field,  describing  the  direction 
of  anisotropy,  and  a  residual  pattern  obtained  by  describing  the  image  in  a  coordinate 
system  built  from  the  flow  field.  The  texture  orientation  was  estimated  from  Laplacian  of 
Gaussian  filters.  Rao  and  Schunck  [62]  proposed  another  algorithm  based  on  the  gradient 
of  the  Gaussian.  Their  new  algorithm  incorporated  a  more  sophisticated  scheme  for 
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computing  the  coherence  of  the  flow  field. 

In  this  chapter,  an  enhancement  algorithm  based  on  multiscale  wavelet  analysis  is 
described.  Features  were  extracted  by  separable  steerable  filters.  A  coherent  image  and 
phase  information  were  then  generated.  A  nonlinear  function,  integrating  coherent  image 
and  phase  information,  was  applied  to  the  transform  coefficients  at  each  level.  An 
enhanced  image  was  obtained  via  an  inverse  wavelet  transform  of  the  modified  coefficients. 
The  novelty  and  advantage  of  this  algorithm  compared  to  existing  techniques  (including 
those  previously  developed  under  this  grant)  lies  in  its  detection  of  directional  features  and 
removal  of  unwanted  perturbations  (artifacts). 

2.4.2  Background 

In  this  section  we  briefly  describe  the  mathematical  background  and  fundamental  ideas 
used  in  subsequent  sections. 


Wavelet  Transforms 

Wavelet  transforms  have  become  well  recognized  as  useful  tools  for  many  applications 
in  signal  processing.  A  function  ip{x)  is  said  to  be  a  wavelet  if  and  only  if  its  Fourier 
transform  'i/>(^)  satisfies  the  admissibility  condition 
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(90) 


This  condition  implies  that 


(91) 


This  means  that  ^^{x)  will  have  at  least  some  oscillations. 

Wavelets  constitute  a  family  of  functions  derived  from  one  single  function  V'  {mother 
wavelet)  by  dilations  and  translations 

(^)  ’ 

where  a  e  R^,b  E  R.  The  idea  of  the  wavelet  transform  is  to  represent  any  function  f{t)  as 
a  superposition  of  wavelets.  The  continuous  wavelet  transform  is  defined  as 

1  /'x  -  b\ 

Wf{a,b)  =  {f,'ipa,b)  =  f{x)i^\^—^yx,  (93) 
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where  il)  denotes  the  complex  conjugate  of  ip.  A  function  can  be  reconstructed  from  its 
wavelet  transform  by  means  of  the  “resolution  of  identity”  formula 
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In  the  practice  of  digital  mammography,  one  prefers  to  express  /  as  a  discrete 
superposition.  Let  us  discretize  the  translation  and  dilation  parameters  of  the  wavelet  in 
Equation  (92): 


=  Oo  "  -  nbo),  (95) 

where  a  =  a^, 6  =  n6oa^,  with  m,neZ,  and  oo  >  1, 5o  ^  0.  On  this  discrete  grid,  the 
wavelet  transform  is  simply 


^  ^“tOO 

Wf{m,n)-aQ^  f{x)ip{aQ'^x-nbo)dx.  (96) 

J  — OO 

The  original  signal  can  be  approximated  as  linear  combinations  of  the  wavelet  bases. 


fix)  ~  Wf{m,n)ipmAx)-  (97) 

mfn 

One  popular  discretization  is  to  choose  Oq  =  2,  bo  =  1, 

'4’m,nix)  =  2"^ ^(2~"^x  -  n),  (98) 


which  are  called  dyadic  wavelets. 


Steerable  Filters 

A  function  f{x,  y)  is  called  “steerable”  if  it  can  be  expressed  as  a  linear  combination  of 
rotated  versions  of  itself.  The  fundamental  idea  of  steerable  filters  is  to  apply  basis  filters 
which  correspond  to  a  fixed  set  of  orientations  and  interpolate  between  each  discrete 
response.  Thus,  one  must  decide  the  number  of  “basis  filters”  and  the  corresponding 
interpolation  functions.  As  defined  in  [20]  a  steering  constraint  may  be  formulated  by 

M 

where  M  is  the  number  of  basis  functions  required  to  steer  a  function  /^'(x,  y). 

Hereafter,  it  will  be  more  convenient  to  work  in  polar  coordinates  r  =  ^Jx^  +  y^  and 
(j)  —  arg(a;,  y).  Let  /  be  any  function  that  can  be  expressed  as  a  Fourier  series  in  polar 
angle  (f)\ 

fir,f))^  an(?’)e^"'^,  (499) 

n^-N 


77 


where  j  = 

The  theorem  below  was  posed  by  [20]  and  is  included  for  the  clarity  of  discussion. 
Theorem  1:  The  steering  condition  (99)  holds  for  a  function  /  expanded  in  the  form  of 
(100)  if  and  only  if  the  interpolations  ki{9)  are  solutions  of 


Then,  f^{r,(j))  may  be  expressed  as 

M 

f{r,  0)  =  ki{0)gi{r,  (j)),  (102) 

i=l 

where  pi(r,  f)  can  be  any  set  of  functions. 

Measure  of  Coherence 

Texture  plays  an  important  role  in  many  machine  vision  and  image  processing  tasks 
including  surface  inspection,  scene  classification,  surface  orientation  and  shape 
determination  [63].  Texture  patterns  may  be  characterized  by  extracting  measurements 
that  quantify  the  nature  and  directions  of  pattern.  Most  breast  carcinomas  have  the 
appearance  of  stellate  lesions  consisting  of  a  central  mass  surrounded  by  radiating  spicules 
[54].  The  spicules  radiate  outward  in  all  directions  and  vary  in  length.  This  provides  an 
important  cue  for  early  cancer  detection. 

Much  attention  has  been  given  to  the  notion  of  decomposing  an  intensity  image  into 
intrinsic  images  to  extract  meaningful  information  [64,  65].  These  intrinsic  properties 
represent  basic  components  of  the  image  formation  process  and  therefore  reveal  features 
“hidden”  inside  an  image.  The  information  they  provide  is  beyond  the  intensity  image 
alone. 

Rao  and  Schunck  [62]  defined  the  orientation  field  of  a  texture  image  to  consist  of  two 
images  —  an  angle  image  and  a  coherence  image.  The  angle  image  denotes  the  dominant 
local  orientation  at  each  point  and  the  coherence  image  represents  the  degree  of  anisotropy 
at  each  point.  They  strongly  advocated  the  use  of  angle  and  coherence  images  as  intrinsic 
images.  In  this  chapter,  we  investigated  the  efficiency  of  these  two  representations  to 
capture  and  enhance  features  of  importance  to  mammography. 

2.4.3  Methodology 

Our  algorithm  consists  of  the  following  four  steps. 
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(1)  Multiscale  Wavelet  Transform:  Wavelet  transforms,  owing  to  their  localization 
characteristics,  are  powerful  tools  of  analysis  for  many  signal  and  image  processing 
applications.  Multiscale  analysis  can  extract  features  at  distinct  scales  and  provide  local 
information  often  hidden  in  an  original  mammogram.  In  our  algorithm,  a  digitized 
mammogram  was  decomposed  using  a  fast  wavelet  transform  algorithm  (FWT)  [7].  In 
order  to  obtain  wavelet  coefficients  at  each  level  without  downsampling,  an 
“algorithme  a  trous"  (algorithm  with  holes)  [13,  66]  was  implemented.  Let  denote  an 
original  image  and  D^f  be  obtained  by  inserting  2*  —  1  zeros  between  every  pair  of  the 
coefficients  representing  /.  {D^f)x  and  {D^f)y  stand  for  carrying  out  convolution 
operations  with  the  filter  D^f  along  x  and  y  directions,  respectively.  The  decomposition 
and  reconstruction  equations  at  level  i  are  as  follows: 

Decomposition: 

=  sW  {D^g)y.  (103) 

Reconstruction: 

S*  =  *  {D^k)^  *  {DH)y  +  *  (DH)^  *  {D^k)y  +  *  (D^h)^  *  {D^h)y,  (104) 

where  indicates  discrete  convolution,  and  h,g,k,  and  I  are  filters  whose  Fourier 
transforms  K(uj),  and  L{u>),  respectively)  satisfy  [7] 


G{uj)K{u;)  +  \H{uj)\^  ^1, 
L{uj)  = - - • 


(105) 


(2)  Separable  Steerable  Filters:  A  filter  is  called  “steerable”  if  the  filter  at  an  arbitrary 
orientation  can  be  expressed  as  a  linear  combination  of  a  set  of  basis  filters,  generated  from 
rotations  of  a  single  kernel  [20].  Steerable  filters  [20,  67,  68],  which  can  be  adaptively 
adjusted  to  arbitrary  orientation,  were  used  to  detect  stellate  patterns  of  spicules  and 
locate  feature  orientations  more  precisely.  As  pointed  out  by  [21],  the  separability  property 
of  the  filters  sped  up  computations  considerably  when  convolved  with  a  large  image  matrix. 
In  our  algorithm,  we  used  three  basis  functions  as  steerable  filters.  The  x-y  separable 
steerable  approximations  of  filter  kernels  were  generated  by  Singular  Value  Decomposition 
(SVD)  [20,  21].  Using  a  set  of  separable  steerable  filters,  local  energy  (M*)  and  associated 
dominant  directions  (A*)  were  determined  by  the  basis  functions  and  their  Hilbert 
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transforms  [20,  69,  70,  71].  We  applied  separable  steerable  filters  to  an  original  image  5° 
and  low-pass  filtered  images  at  each  level  5*,  i  =  1, . . .  n,  to  capture  salient  multiscale 
features.  Discrete  realizations  of  the  three  basis  functions  applied  in  our  algorithm  are 
shown  in  Figure  28. 


Figure  28:  Two-dimensional  impulse  responses  of  three  basis  functions. 


(3)  Coherence  Maps:  A  coherence  map  is  an  image  showing  a  local  measure  of  the 
degree  of  anisotropy  of  flow  [61,  62].  If  the  orientations  of  a  texture  pattern  at  any  point 
{xi,yi)  are  coherent,  then  magnitude  and  phase  information  are  important  and  should  be 
emphasized.  Conversely,  if  the  orientations  are  not  coherent,  the  magnitude  and  phase 
information  can  be  neglected  or  attenuated.  Kass  and  Wilkin  [61]  suggested  a  simple  way 
of  measuring  strength  of  coherence  by  finding  the  ratio 


X{h  k) 


(106) 


where  W {j,  k)  denotes  a  local  weighting  function  with  unit  integral,  J {j,  k)  denotes  the 
squared  gradient  vector  at  {j,k),  and  j  •  |  denotes  absolute  value. 

An  alternative  measure  of  coherence  was  proposed  by  Rao  and  Schunck  [62]  and  was 
obtained  by  weighting  the  energy  with  the  normalized  projection  of  energy  within  a 
specified  window  (W)  onto  the  central  point  {j,  k)  of  the  window.  The  coherence  (C*)  was 
expressed  as 


C%k)  =  M\j,k) 


|M*(m,  n)cos(A*(j,  A;)  —  A*(m,n))] 

(m,n)GW 

(m,7i)€W 


(107) 


where  M^{j,k)  and  A^{j,k)  denote  energy  and  phase  of  point  {j,k)  at  level  i,  respectively. 
This  coherence  measure  incorporated  the  gradient  magnitude  and  hence  placed  more 


80 


Level  0 


Level  I 


Level  2 


Figure  29:  Overview  of  processing  for  Steps  1-3. 

weight  on  regions  that  had  higher  visual  contrast.  This  method  performed  better  than 
previous  measures  applied  to  the  similar  data  [62].  In  our  algorithm,  we  implemented  this 
method  to  measure  the  coherence  map  at  each  level.  A  5  x  5  window  was  used  to  perform 
smoothing.  The  measure  of  coherence  C*  was  obtained  from  Equation  (107).  The 
coherence  and  orientation  data  extracted  the  more  salient  features  of  spiculated  lesions.  A 
schematic  diagram  of  Steps  1-3  is  shown  in  Figure  29. 

(4)  Nonlinear  Operators:  So  far,  we  have  computed  all  the  information  we  need  in  our 
algorithm.  A  nonlinear  operation  was  then  applied  within  each  level  to  precisely  modify 
transform  coefficients.  This  operation  integrated  both  coherence  map  and  phase 
information. 

A.  Modification  from  Coherence  Map:  Let  C^{j,k)  denote  the  coherence  measure  of  point 
{j,  k)  at  some  level  i.  Modifications  from  coherence  measure  were  obtained  by  a  nonlinear 
function  expressed  as 

MC\j,k)  =  ,if  C\j,k)^0, 

=  0  ,if  C%k)=0, 

such  that  a  coefficient  was  emphasized  if  its  coherence  measure  was  large  and  attenuated  if 
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small. 

B.  Modification  from  Phase:  Phase  information  is  important  to  characterize  well  an 
oriented  texture  and  therefore  we  did  not  neglect  its  contribution  to  the  modification  of 
coefficients.  We  applied  a  sinusoidal  weighting  to  the  phase  information.  The  “detail” 
sub-bands  of  wavelet  coefficients  obtained  in  Step  1  included  two  components:  the 
component  along  x  direction  and  the  component  along  y  direction.  The  x  component  was 
obtained  by  high-pass  filtering  along  x  direction,  hence  mostly  vertical  features  within  the 
mammogram  were  detected.  We  emphasized  the  points  whose  dominant  orientations  were 
near  0  and  tt  (with  respect  to  the  vertical  axis).  Thus,  the  modification  from  phase 
information  was 


MA^{j,  k)  =  0.01  +  I  cos(^*(i,  k))\.  (108) 

Note  that  a  constant  (0.01)  was  added  to  the  above  modification  so  that  the  phase  factor 
would  not  be  neglected  when  the  phase  was  equal  to  |  or 

The  y  component  was  obtained  by  high-pass  filtering  along  the  y  direction,  hence 
horizontal  features  of  an  image  were  detected.  We  emphasized  the  points  whose  dominant 
orientations  were  near  |  and  ^  (with  respect  to  the  vertical  axis).  The  modification  from 
phase  information  was 

MA\j,  k)  =  0.01  +  I  sin(W(j,  k))\.  (109) 

Note  that  a  constant  (0.01)  was  added  to  the  above  modification  so  that  the  phase  factor 
would  not  be  neglected  when  the  phase  was  equal  to  0  or  tt. 

The  modifications  from  coherence  map  and  phase  at  level  i  —  1  were  combined  to 
adjust  the  wavelet  coefficients  at  level  i.  The  final  modification  was  therefore 

MW*  =  •  W^  (110) 

where  T*  was  a  constant  at  each  level.  A  schematic  diagram  of  the  nonlinear  operator  at 
level  i  is  shown  in  Figure  30.  These  modified  coefficients  were  then  reconstructed,  via  an 
inverse  fast  wavelet  transform,  to  enhance  the  visualization  of  possible  lesions. 

2.4.4  Experimental  Results 

Our  algorithm  was  applied  to  enhance  mammograms  using  oriented  information.  In  order 
to  capture  distinct  directions  of  subtle  features,  steerable  filters  were  used.  Figure  31  shows 
the  capability  of  steerable  filters  for  detecting  and  enhancing  objects  of  distinct  direction 
and  scale.  The  phantom  image  contained  24  diamond  shaped  objects  with  four  directions 
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Level  i 


Figure  30:  A  schematic  diagram  of  nonlinear  operator  at  level  i. 

(0°,30°,60°,  and  90°)  and  six  scales,  as  shown  in  Figure  31(a).  After  applying  our 
algorithm  to  this  image,  the  borders  of  these  objects  at  each  scale  were  clearly  enhanced. 
The  phantom  objects  in  the  last  two  rows  of  the  image,  were  very  close  to 
microcalcifications.  This  demonstrated  the  ability  of  this  algorithm  to  detect 
microcalcification  clusters  in  mammograms,  without  artifacts. 

Three  examples  of  malignant  lesions  with  distinct  radiographic  signs  of  cancer  were 
processed  to  show  the  effectiveness  of  our  algorithm.  The  images  were  of  matrix  size 
(512  X  512).  For  each  case,  both  global  and  regions  of  interest  (ROI)  were  shown  along 
with  the  corresponding  enhanced  ROI  image. 

Calcifications 

Figure  32  shows  a  mammogram  (mamOOQlml)  with  microcalcification  clusters.  The 
original  mammogram  is  shown  in  Figure  32(a).  Figure  32(b)  shows  an  original  suspicious 
area  of  the  mammogram.  After  enhancement,  clusters  of  calcifications  appear  clearly  in  the 
center  of  the  image. 

Stellate  lesions 

A  mammogram  (mam041rcc)  with  a  stellate  lesion  is  shown  in  Figure  33.  Figure  33(a) 
shows  the  original  image.  Figure  33(b)  presents  an  original  digital  radiograph  with  a 
partially  obscured  irregular  mass  in  the  center  of  the  image  matrix.  After  applying  our 
algorithm  to  this  image,  the  enhanced  image  makes  obvious  spiculated  lesions  around  the 
mass.  The  mass  itself  was  also  enhanced,  as  shown  in  Figure  33(c).  These  clear  spiculated 
patterns  suggested  radiologists  that  this  mass  is  more  likely  malignant,  rather  than  benign. 
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Figure  31:  The  ability  of  detecting  objects  oriented  along  distinct  directions  at  various  scales: 
(a)  Original  test  image  (mathematical  phantom) ,  (b)  enhanced  image  after  processing. 

Masses 

A  mammogram  (mam0041cc)  with  a  mass  tumor  is  shown  in  Figure  34.  The 
craniocaudal  view  of  the  left  breast  shown  in  Figure  34(b)  shows  an  irregular  spiculated 
mass  in  retroglandular  fat.  The  enhanced  version  shown  in  Figure  34(c)  better  delineates 
the  margins  of  the  mass. 

2.4.5  Summary 

An  enhancement  algorithm  relying  on  multiscale  wavelet  analysis  and  extracting  oriented 
information  at  each  scale  of  analysis  was  investigated.  The  evolution  of  wavelet  coefficients 
across  scales  characterized  the  local  shape  of  irregular  structures.  Using  oriented 
information  to  detect  the  features  of  an  image  appears  to  be  a  promising  approach  for 
enhancing  complex  and  subtle  structures  of  the  breast.  Steerable  filters  which  can  be 
rotated  at  arbitrary  orientations  can  reliably  find  visual  cues  within  each  spatial-frequency 
sub-band  of  an  image.  “Coherence  measure”  and  “dominant  orientation”  clearly  helped  us 
to  discriminate  features  from  complex  surrounding  tissue  typical  in  mammograms. 

Existing  and  previous  multiscale  enhancement  approaches  [72,  32,  59]  attempted  to 
enhance  an  image  by  detecting  edges.  Unfortunately,  most  edge  detection  algorithms  can 
not  distinguish  between  “authentic”  edges  and  phantom  edges.  In  contrast,  this  algorithm 
relied  upon  a  coherence  map  and  phase  information  which  resulted  in  an  enhancement 


(b)  (c) 


Figure  34:  Mammogram  with  a  mass:  (a)  Original  mammogram,  (b)  ROI  image,  (c)  en¬ 
hanced  ROI  image. 
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naturally  close  to  the  original  image.  This  type  of  artifact  free  enhanced  image  can  provide 
more  “obvious”  (and  familiar)  visual  cues  for  radiologists. 

2.5  Multivoice  Undecimated  Wavelet  Transforms 

2.5.1  Introduction 

For  many  analysis/synthesis  applications  the  sampling  provided  by  the  traditional  wavelet 
transform  is  inadequate.  In  particular,  a  finer  grid  is  needed  when  computing  the  wavelet 
transform  for  applications  that  require  “good”  time-frequency  localization.  This  can  be 
achieved  by  computing  an  undecimated  wavelet  transform  and  the  addition  of  voices.  In 
this  chapter  we  review  the  tools  that  will  allow  us  to  compute  an  exact  multivoice 
undecimated  wavelet  transform  for  an  arbitrary  wavelet.  Furthermore,  we  introduce  a 
wavelet  function  that  exhibits  nearly  optimum  time-frequency  localization,  a  desirable 
property  in  signal  analysis  as  we  will  see  below. 

The  chapter  is  organized  as  follows:  In  Sections  2.5.2  and  2.5.3  we  review  the 
short-time  Fourier  transform  and  the  wavelet  transform,  including  the  theory  of  frames, 
which  will  allow  us  to  carry  out  wavelet  analysis/synthesis  for  an  arbitrary  wavelet.  In 
Sections  2.5.4  and  2.5.5  we  review  the  a  trous  and  Mallat’s  algorithms  and  learn  that 
although  the  a  trous  algorithm  was  originally  devised  as  a  computationally  efficient 
implementation  of  the  wavelet  transform,  it  is  more  properly  viewed  as  a  nonorthogonal 
multiresolution  decomposition  for  which  the  discrete  wavelet  transform  computes  a 
sampled  continuous  wavelet  transform  exactly.  In  Sections  2.5.6  and  2.5.7  we  introduce  the 
sine-Gabor  wavelet  and  show  that  it  exhibits  nearly  optimum  time-frequency  localization 
and  satisfies  the  conditions  of  a  wavelet  frame.  Finally,  in  Section  2.5.8,  we  summarize  our 
results  and  describe  related  ongoing  work. 

2.5.2  The  Short-Time  Fourier  Transform 

The  goal  of  signal  analysis  is  to  extract  relevant  information  from  a  signal  by  transforming 
it  into  a  representation  where  the  properties  of  the  signal  are  more  evident.  For  the 
analysis  of  stationary  signals,  that  is,  signals  whose  properties  do  not  evolve  with  time,  one 
example  of  such  representation  is  the  Fourier  transform 

/OO 

dt. 

'OO 

The  Fourier  transform  can  be  viewed  as  the  inner  product  of  the  signal  s{t)  and  the 
sinusoidal  wave  The  analysis  coefficient  s{uj)  -  {J^s){u)  measures  the  “strength”  of 
the  sinusoidal  wave  of  frequency  uj  in  the  signal  s{t). 


88 


Fourier  analysis  works  well  for  signals  composed  of  a  few  stationary  components. 
However,  any  abrupt  change  in  time  in  a  non-stationary  signal  s{t)  is  spread  out  over  the 
whole  frequency  axis  in  s(a;).  In  order  to  adapt  the  Fourier  transform  to  non-stationary 
signals,  Gabor  [73]  introduced  a  new  transform  by  using  a  windowed  function  w{t)  in  the 
Fourier  integral 

/OO  _ 

s{t)w{t  —  dt. 


A  function  w{t)  qualifies  as  a  window  function  if  it  is  possible  to  identify  its  center  and 
standard  deviation  (or  root  mean  square  (RMS)  duration)  defined  by 


and 


m{w)  = 


i|i(;(t)P  dt 


1  (  r°°  1 

^  f^XJ  ,  (111) 

respectively  [74].  In  the  original  Gabor  transform,  the  window  function  was  a  Gaussian 
[73],  however,  the  transform  is  valid  for  any  type  of  window  function  and  is  traditionally 
referred  to  as  the  short-time  or  window  Fourier  transform. 

The  short-time  Fourier  transform  can  be  viewed  as  the  inner  product  of  the  signal  s{t) 
with  the  family  of  functions 


WT,u}{t)  =  w{t  —  G 


that  is. 


(^^s)(r,a;)  =  (s(t),Wr,u;(t)} .  (112) 

It  is  easy  to  show  that  Wr,ui{t)  is  a  window  function  with  center  and  standard  deviation 
given  by 


and 


=  m{w)  +  T 

(113) 

cr{Wr,c)  =  cr{w), 

(114) 

respectively.  From  Equations  112,  113  and  114  we  see  that  in  the  time  domain  the  analysis 
coefficient  (^«,s)(r,a;)  essentially  depends  on  the  values  of  s{t)  for 

t  G  [miw)  -f-  r  —  (y{w),  m{w)  t  +  o{w)\ .  This  is  called  time  localization.  It  follows  from 
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the  above  observation  that  two  pulses  in  time  can  be  discriminated  only  if  they  are  more 
than  2a{w)  apart.  This  is  referred  to  as  the  time  resolution  of  the  short-time  Fourier 
transform. 

By  applying  Parseval’s  theorem  to  Equation  112  we  obtain  an  alternative  formula  for 
the  short-time  Fourier  transform 

{QwS){t,U))  =  ^  (115) 

where 


^r,uj(7)  =  w(7  -  {t,uj)  e  R^. 

Suppose  w{'y)  =  is  also  a  window  function  with  center,  m{w),  and  standard 

deviation,  a{w).  Then  it  is  easy  to  show  that  Wr,uj{7)  is  a  window  function  with  center  and 
standard  deviation  given  by 

'm{wT,uj)  =  fn{w)  +  LO  (116) 

and 

o-(u)t,w)  =  cr(u)),  (117) 

respectively.  From  Equations  115,  116  and  117  we  see  that  in  the  frequency  domain  the 
analysis  coefficient  (^^s)(r,  cu)  essentially  depends  on  the  values  of  s('y)  for 
7  G  [m(w)  +(jj  -  (7(w),  m(w)  T-u)  +  a{w)] .  This  is  called  frequency  localization.  It  follows 
from  the  above  observation  that  two  pure  sinusoids  can  be  discriminated  only  if  their 
frequencies  are  more  than  2a (w)  apart.  This  is  referred  to  as  the  frequency  resolution  of 
the  short-time  Fourier  transform. 

From  Equations  112  through  117  we  observe  that  {Gws){t,uj)  yields  a  time-frequency 
representation  of  s{t).  The  analysis  coefficient  {Gws){t,u)  depends  on  the  time-frequency 
window 

[m{w)  +  T  -  a{w),m{w)  +  t  -\-  a{w)]  x  [m{w)  -I-  a;  -  a{w),m{w)  +  co  +  cr(ry)] . 

Increasingly  accurate  localization  in  time  and  frequency  is  not  possible  because  the  area  of 
the  time-frequency  window  is  lower  bounded  by  2,  that  is 

4a{'w)a{w)  >  2.  (118) 

This  is  referred  to  as  the  uncertainty  principle  or  Heisenberg  inequality.  Equality  is 
satisfied  if  and  only  if  the  window  function  is  a  Gaussian  as  in  the  Gabor  transform. 
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The  short-time  Fourier  transform  is  an  isometry  (to  a  proportionality  coefficient)  from 
L2(E)  into  L2(E2)  [75],  that  is 

1  p+oo  p+oo 

The  function  s{t)  is  reconstructed  from  {GwS){t,lo)  with  the  formula 

-j  ^+oo  p-\-oo 

=  o  II  112  /  /  (Mir, ^)w{t  -  dTduj. 

27r\\w\\l  y_oo  7-00 

The  short-time  Fourier  transform  is  a  redundant  representation.  Instead  of  computing 
{Gws){r,u)  for  all  values  {r,u))  in  it  is  possible  to  uniformly  sample  both  r  and  w  such 
that  the  representation  is  complete  and  stable  [75].  Let  tq  and  loq  be  the  sampling  intervals 
in  the  time  and  frequency  domains  respectively.  Then,  the  discrete  short-time  Fourier 
transform  is  defined  by  [75] 

(G^s)^,n  =  {Qws){mTo,n(Jo),  (m,n)  G 

To  reconstruct  any  function  s{t)  E  L‘^{R)  from  the  set  of  samples  (m,n)  G 

the  operator 

L2(R)  ^ 

must  be  invertible  on  its  range  and  have  a  bounded  inverse  [75].  In  order  to  invert  G^, 
Daubechies  [76]  has  shown  that  tq  and  uq  must  verify 

Cc^o'^O  ^  ^TT. 

A  drawback  of  the  short-time  Fourier  transform  is  that  once  a  window  function  has 
been  chosen  the  time-frequency  resolution  of  the  transform  is  fixed  over  the  entire 
time-frequency  plane.  It  is  therefore  impossible  with  the  short-time  Fourier  transform  to 
analyze  at  the  same  time  transient  signal  components  with  good  time  resolution  and 
quasi-stationary  signal  components  with  good  frequency  resolution  [77]. 


2.5.3  The  Wavelet  Transform 


The  wavelet  transform  overcomes  the  resolution  limitation  of  the  short-time  Fourier 
transform  by  letting  the  time  and  frequency  resolution  vary  in  the  time-frequency  plane. 
The  wavelet  transform  is  defined  by  [74] 


{W^s){a,b)  =  |o|“^  J  s{t)il)  ^  o.,bER,  {a  ^  0).  (119) 
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In  order  to  reconstruct  s{t)  from  its  wavelet  transform,  the  Fourier  transform  of  must 
satisfy  [74] 


0,1 


/+00 

■oo 


l'0(7)P 

l7l 


dj  <  oo. 


(120) 


A  function  is  said  to  be  a  basic  wavelet  if  it  is  possible  to  reconstruct  s{t)  from  its 
wavelet  transform. 

For  any  satisfying  Equation  120,  the  wavelet  transform  is  an  isometry  (to  a 
proportionality  coeflBcient)  from  L^(E)  into  L^(E^)  [74],  that  is 


12  _  1 

0, 


|(>V^s)(a,6)|^  ^db. 


J  — OO  J  — OO 

The  function  s{t)  is  reconstructed  from  (>V^s)(o,  6),  a,  6  G  E,  (a  7^  0),  with  the  formula  [74] 


^  /r  /-i” (^) } 


^db. 


The  wavelet  transform  can  be  viewed  as  the  inner  product  of  the  signal  s{t)  with  the 
family  of  functions 


1  l-i  / 

ft-b\ 

|a|  2-0 

a 

,  a,  6  G  E,  (a  /  0), 


(121) 


that  is, 


{yV,ps){a,b)  =  {s{t),'ipa-,b{t))  ■ 


(122) 


Suppose  is  a  window  function  with  center  m{'ip)  and  standard  deviation  a{'tp). 
Then,  it  is  easy  to  show  that  '4>a-,b{t)  is  a  window  function  with  center  and  standard 
deviation  given  by 


'm{'ipa-,b)  =  am{ij;)  +  b  (123) 

and 

cr{.i^a-,b)  =  i«k(^),  (124) 

respectively.  From  Equations  122,  123  and  124  we  have  that  in  the  time  domain  the 
analysis  coefficient  (>V^s)(a,  6)  essentially  depends  on  the  values  of  s{t)  for 
t  G  [am{il>)  +  b-  \a\(j{ijj),  am{'ip)  +  b+  |a|(T(^)] .  It  can  be  verified  [74]  that  if  a  basic 
wavelet,  is  also  a  window  function  then  is  necessarily  in  L^(E)  so  that  its  Fourier 
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transform,  '0(w)  is  a  continuous  function.  It  follows  from  Equation  120  that  '0(u’)  must 
vanish  at  the  origin,  that  is, 


dt  =  0. 


(125) 


By  applying  Parseval’s  theorem  to  Equation  122  we  obtain  an  alternative  formula  for 
the  wavelet  transform 


(126) 


where 


0,6  G  E,  {a^  0). 


Suppose  '0(t)  also  a  window  function  with  center,  Tniijj),  and  standard  deviation,  (7('0). 
Then  it  is  easy  to  show  that  '>pa-,b{'y)  is  a  window  function  with  center  and  standard 
deviation  given  by 


(127) 


and 


a{^p) 


O 


(128) 


respectively.  From  Equations  126,  127  and  128  we  have  that  in  the  frequency  domain  the 
analysis  coefficient  (VV-0s)(a,  6)  essentially  depends  on  the  values  015(7)  for 

m{'ip)  cr(Vj)  m('ip)  

a  |a|  ’  a  \a\ 

From  Equations  122  through  128,  {W^s){a,b)  yields  a  time-frequency  representation 
of  s{t).  The  analysis  coefficient  {W^s){a,b)  depends  on  the  time-frequency  window 


je 


[am{'ip)  +  b—  \a\a{'ip),  +  b+  |a|(r('i/’)]  x 


m('0)  <j('0)  m('0)  _  (j('0) 

a  lol  ’  o  |o| 


In  contrast  to  the  short-time  Fourier  transform,  both  the  time  and  frequency  resolution  in 
the  wavelet  transform  vary  in  the  time-frequency  plane.  The  following  properties  derive 
from  the  time-frequency  localization  characteristics  of  the  wavelet  transform: 


1.  The  ratio  of  the  standard  deviation,  a{'^)/\a\,  of  the  frequency  window  and  its  center 
frequency,  m('ip)/a,  is  given  by  which  is  independent  of  the  location  of  the 

myyj) 

center  frequency.  This  is  called  constant-Q  frequency  analysis. 
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2.  The  time-frequency  window  narrows  in  time  and  widens  in  frequency  for  large  center 
frequency  m{-ip)/a  (small  a  >  0),  and  widens  in  time  and  narrows  in  frequency  for 
small  center  frequency  m{'ip)la  (large  a  >  0). 


Recall  that  the  wavelet  transform  in  Equation  119  can  be  viewed  as  the  inner  product 
of  s{t)  with  the  family  of  functions  ipa-,b{i)  defined  in  Equation  121.  Notice  that  this  family 
of  functions  is  obtained  from  a  single  basic  wavelet  by  translating  it  by  h  and  dilating  it  by 
a.  The  basic  wavelet  is  usually  referred  to  as  the  mother  wavelet.  The  constant  |a|  2  in 
Equation  121  is  used  for  energy  normalization  purposes,  that  is  \\'ipa-,b{t)\\2  —  Il'0(^)ll2- 
general  the  mother  wavelet  'ip{t)  is  normalized  to  have  its  energy  equal  to  one,  that  is 
||-0(t)||2  =  1.  We  assume  this  restriction  from  now  on. 

In  practice  the  parameter  a  is  restricted  to  a  >  0,  so  that  the  center  frequency  of 
'^a\b{l)  in  Equation  127  is  restricted  to  positive  frequencies.  In  this  context  the  parameter 
a  is  usually  referred  to  as  the  scale  parameter.  In  order  to  reconstruct  s{t)  from  its  wavelet 
transform  restricted  to  a  >  0,  the  Fourier  transform  of  •0(t)  must  satisfy  a  more  strict 
admissibility  condition  than  Equation  120  given  by  [74] 

|?^(7)|2  _  _ 

■ .  "  ■  2' 


7 


■  dy  = 


r'+oo  1 

■dj=-C^<oo. 


7 


(129) 


For  any  V’(^)  satisfying  Equation  129  we  have  the  following  reconstruction  formula  [74] 


+00  /»  +  oo 


0-0  J  —QQ  J  0 


{W^s){a,b) 


The  wavelet  transform  is  a  redundant  representation.  Instead  of  computing 
(>V^s)(a,  b)  for  6  G  M  and  a  G  E+,  it  is  possible  to  choose  an  exponential  sampling  of  the 
scale  parameter,  a  =  a^,  j  G  Z,  oq  ^  0,  such  that  the  representation  is  complete  and  stable. 
Of  particular  interest  is  the  sequence  of  scales  where  the  elementary  dilation  step  oq  =  2 
because  it  leads  to  an  octave  by  octave  partitioning  of  the  frequency  domain. 

In  order  to  reconstruct  s{t)  from  its  wavelet  transform  restricted  to  6  G  M  and 


-  J 


ai,  j  eZ,  the  Fourier  transform  of  -ip{t)  must  satisfy  a  more  strict  condition  than 


Equation  129  given  by  [74] 


+00 


<  B, 


(130) 


j=-oc 


for  almost  all  7  G  E  for  some  constants  A  and  B  with  0  <  A  <  B  <  00.  Reconstruction  is 
possible  by  using  the  following  formula  [74] 

+00 


j=-oo 


t  —  b 


db, 
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where  the  Fourier  transform  of  is  given  by 


j=-oo 

If  ■ip{t)  satisfies  the  so-called  “stability  condition”  in  Equation  130  with  ao  =  2,  then  ^(t)  is 
called  a  dyadic  wavelet  [74]. 

In  addition  to  sampling  the  scale  parameter  a  =  a^,  j  e  Z,  it  is  possible  to  choose  a 
sampling  of  the  translation  parameter  b  and  still  obtain  a  complete  and  stable 
representation.  Intuitively,  the  discretization  of  b  has  to  be  chosen  such  that  ipaibit) 

“covers”  the  entire  time  axis  for  each  scale  a  =  ai,  j  G  Z.  From  Equation  124  we  see  that 
the  standard  deviation  of  ipa-,b{t)  is  [ol  times  the  standard  deviation  of  Hence,  in  order 
for  'ipa-,b{t)  to  “cover”  the  entire  time  axis  at  scale  a  =  a^,  the  discretization  of  b  has  to  be 
proportional  to  |a|.  It  follows  that  at  scale  a  =  al  the  translation  has  to  be  6  =  kboal, 
where  bo  >  0  is  the  elementary  translation  step.  This  leads  to  a  discretized  family  of 
wavelets  given  by 

=  O'o^'tp  ^ 

In  order  to  reconstruct  s(t)  from  its  wavelet  transform  restricted  to  6  =  kboa^,  j,k  eZ 
and  a  =  al,  j  e  Z,  must  satisfy  a  stricter  condition  than  Equation  130.  The 
admissibility  condition  for  this  reconstruction  is  the  existence  of  H  >  0,  H  <  oo  so  that 

^Ikll2  <  <  B\\s\\l 

j,k^Z 

for  all  s{t)  G  I/^(R)  [8].  In  other  words,  the  family  of  functions  defined  in  Equation 

131  constitute  a  frame.  A  brief  review  of  frames  is  presented  at  the  end  of  this  section. 

Given  a  frame  of  wavelets,  reconstruction  of  s{t)  from  the  {s{t),'ipj^k{t))  is  possible  by 
using  the  following  formula 

+  CX)  H-oo 

^W=E  (132) 

oo  k=—oo 

where  the  family  of  functions  j,k  E  Z  is  called  the  dual  frame  of  j,k  E  Z 

(see  the  section  below). 

If  the  dual  frame  of  j,k  E  Z  is  of  the  form 

j,k  E  Z, 
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where 


^o,k{t)  =  ao  j  ^  -  ^^o) , 

then  Equation  132  is  referred  to  as  a  wavelet  series  and  'ip^t)  is  known  as  the  dual  wavelet 
of  ?/’(t).  Given  a  wavelet  frame,  the  existence  of  a  dual  wavelet  is  important  for 
computational  purposes.  Below,  we  overview  the  classification  of  wavelets  according  to 
orthogonality  and  highlight  the  existence  of  a  dual  wavelet  in  each  case.  A  wavelet  can  be 
classified  according  to  orthogonality  as  follows: 

1.  A  wavelet  'ip{t)  is  called  semiorthogonal  if  it  satisfies  =  0, 

j  7^  ^5  j:  k,l,7n  G  Z. 

Remarks:  It  can  be  shown  [74]  that  for  every  semiorthogonal  wavelet,  ip{t),  there 
exist  a  dual  wavelet,  such  that  the  pair  satisfies  the  biorthogonality 

property 

=  Sj^iSk,mJ,k,l,m  G  Z.  (133) 

2.  A  wavelet  ip{t)  is  called  nonorthogonal  if  it  is  not  a  semiorthogonal  wavelet. 

3.  A  semiorthogonal  wavelet  'ip{t)  is  called  orthogonal  if  it  satisfies 

jjk, 1,171  G  Z. 

Remarks:  It  can  be  shown  that  an  orthogonal  wavelet  is  self  dual. 

4.  A  nonorthogonal  wavelet  ip^t)  is  called  biorthogonal  if  there  exists  a  dual  wavelet, 

ip{t),  such  that  the  pair  satisfies  the  biorthogonality  property  in  Equation 

133. 

The  computation  of  Equation  132  can  be  a  burden  for  nonorthogonal  wavelets  for 
which  a  dual  wavelet  does  not  exist.  In  such  a  case,  it  is  advantageous  to  work  with  frames 
which  are  almost  tight  (“snug  frames”),  i.e.,  frames  which  have  B/A  -  1  <C  1,  because  the 
can  be  approximated  by  'ipj,k{t)  (see  the  section  below  and  [8]  for  more  details). 

Frames 

In  this  section  we  present  a  brief  review  of  frames.  For  a  more  detailed  and  rigorous 
account  of  frames  see  [8]. 
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A  family  of  functions  (pj,  j  €  J  in  a  Hilbert  space  H  is  called  a  frame  if  there  exist 
A  >  0,  H  <  oo  so  that  for  all  f 

A\\fr<Y^\{f,V,)f<B\\f\\\  (134) 

where  A  and  B  are  called  the  frame  bounds.  If  the  two  frame  bounds  are  equal,  the  frame 
is  called  a  tight  frame.  In  a  tight  frame  we  have,  for  all  f  eU, 

^\u,<Pi)f  =  Awnf, 

jej 


which  implies 


jeJ 


or  (at  least  in  the  weak  sense  [8]), 

f  =  (W5) 

ieJ 

Although  Equation  135  is  reminiscent  of  the  expansion  of  /  into  an  orthonormal  basis,  it  is 
important  to  note  that  frames,  even  tight  frames,  are  not  (orthonormal)  bases.  The  family 
of  functions  tpj,  j  e  S  are  typically  not  linearly  independent.  If  the  frame  is  tight,  and  if 
=  1  for  all  j  G  J,  then  A  =  B  gives  the  “redundancy  ratio.”  If  this  ratio  equals  to  1, 
then  the  tight  frame  is  an  orthonormal  basis  [8]. 

Equation  135  gives  a  trivial  way  to  recover  /  from  the  (/,  <Pj),  if  the  frame  is  tight. 
Consider  now  recovering  /  from  frames  that  are  not  tight.  Let  us  define  the  frame  operator 
F  from  V.  to  as 


Since  (pj  constitute  a  frame,  it  follows  from  Equation  134  that  HE/p  <  5||/|p,  that  is,  F 
is  bounded,  which  means  it  is  possible  to  find  its  adjoint  operator  F*  [8].  The  adjoint 
operator  F*  of  F  can  be  computed  from  the  following  relation 

{F‘cJ)  =  {c,Ff)  =  '£,c,ifPP~) 

ieJ 

=  'y  ^  Cj  {fPji  f) , 
jel 


so  that 


F*c  =  Cj(pj, 
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at  least  in  the  weak  sense  [8].  From  this  it  follows  that 

E  </'  Vi)  V,  =  F’Ff. 

jeJ 

Thus,  the  definition  of  F  implies 

J2^{f.n)f‘=\\FJf={F’Ff,f). 

jei 

Hence,  in  terms  of  F,  the  frame  condition  in  Equation  134  can  be  written  as 

AI  <  F*F  <  BI,  (136) 

where  I  is  the  identity  operator  [8].  This  implies,  in  particular,  that  F*F  is  invertible  (see 
Lemma  3.2.2  in  [8]),  and  that  the  operator  {F*F)~^  satisfies, 

B-^I  <  {F*F)-^  <  A-^L 

Applying  the  operator  {F*F)~^  to  the  family  of  functions  ipj,  j  G  J,  leads  to  another 
family  of  functions  denoted  by  (pj,  j  G  J  where 

p,  =  (F*F)-V„  j  e  J- 

The  family  of  functions  pj,  j  G  J  also  constitutes  a  frame  with  frame  bounds  B~^  and 
[8],  that  is 

It  can  be  verified  (see  Proposition  3.2.3  in  [8])  that  the  associated  frame  operator  F  from 
n  to  ^2(j),  {Ff)j  =  {f,pj)  satisfies  F  =  F{F*F)-\  F*F  =  {F*F)-\  F*F  =  I  =  F*F  and 
PF*  =1  FF*  is  the  orthogonal  projection  operator,  in  F(J),  onto  Ran(F)  =  Ran(F).  The 
family  of  functions  pj,  j  e  3  is  called  the  dual  frame  of  pj,  j  G  J.  It  is  easy  to  verify  that 
the  dual  frame  of  pj,  j  E  3  is  pj,  j  G  J.  From  F*F  =  I  =  F*F  it  follows  that 

=  f  (137) 

jGJI  JGJ 

Hence,  we  have  a  way  to  recover  /  from  the  (/,  Pj),  where  the  only  thing  we  need  to  do  is 
to  compute  the  Pj  =  {F*F)~'^pj.  Note  that,  in  general,  if  the  frame  is  redundant,  there 
exist  other  functions  in  Ti  that  could  equally  well  play  the  role  of  the  pj  and  lead  to  a 
reconstruction  formula.  This  follows  from  the  fact  that  the  pj  are  not  linearly  independent 
in  the  general  case.  Equation  137,  however,  yields  the  most  “economical”  representation  of 
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/  in  the  following  sense  (see  Proposition  3.2.4  in  [8]):  Consider  the  first  half  of  Equation 
137,  /  =  Ejej  (/>  Vj)  H  /  =  Eiej  (/.  some  other  family  of  functions  j  €  J, 

then  it  can  be  shown  [8]  that  E^gji  ^  Similarly,  consider  the 

second  half  of  Equation  137,  /  =  E^eJ  (/’  ^  ^ 

if  not  all  Cj  equal  (/,  <^j),  then  it  can  be  shown  [8]  that  EjGJ  -  Sjej  I I  • 
Computing  the  ipj  involves  the  inversion  of  F*F.  If  B  is  close  to  A 
(r  =  5/y4  -  1  <  1),  then  from  Equation  136  we  have  that  F*F  is  “close”  to  ^  This 
implies  that  {F*F)^^  is  “close”  to  and  that  ipj  is  “close”  to  [8].  This 

motivates  the  following  reconstruction  formula  for  /  [8], 


where 


R  =  I-  --^F*F.  (139) 

It  follows  from  Equations  136  and  139  that  -|^  I  <R<^I-  This  implies  that 
||7^||  <  1^  =  ^  <  1.  If  r  is  small,  the  rest  term  Rf  in  Equation  138  can  be  dropped, 
leading  to  a  reconstruction  formula  for  /  which  is  accurate  up  to  an  L^-error  of  2^11/11  [3]- 
Even  if  r  is  not  small,  it  is  possible  to  write  an  algorithm  for  the  reconstruction  of  /  with 
exponential  convergence  [8].  From  Equation  139  we  see  that 

F’F=^^{I-R). 


This  implies  that 

(F-F)-'  =  j^{I-R)-K 

Since  ||i?||  <  1,  the  series  E^o-^^  converges  in  norm,  and  its  limit  is  (/  -  R)-^  [8].  It 
follows  that 

2  °° 
k=0 

Note  that  the  zeroth  order  term  in  the  above  equation  leads  exactly  to  Equation  138  with 
the  rest  term  dropped.  Better  approximations  can  be  obtained  by  truncating  after  N  terms 
[8], 


2 

A  +  B 


N 
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with 


f  -^{f,  <fj) 

jeJ 


< 


ll/I 


which  becomes  exponentially  small  as  N  increases,  since  ^  <  1.  In  particular,  the  (fj  can 
be  computed  by  an  iterative  algorithm  [8], 

Similarly,  /  can  be  computed  by  an  iterative  algorithm  [8], 


with 


f={F*F)-\F*F)f=\im  h, 

N^oo 


In  —  /jv-i  +  ^  [(/’  ‘^i)  ~  {In-i,  ^j)]  ^j- 


Wavelet  Frames 

Not  all  choices  for  '0(t),  oq,  bo  lead  to  frames  of  wavelets,  even  if  'ip{t)  is  admissible.  In 
[8],  Daubechies  gives  some  general  conditions  on  uq,  bo  under  which  a  frame  is 
obtained  and  derives  estimates  for  the  frame  bounds.  These  results  can  be  summarized  as 
follows  (see  Proposition  3.3.2  in  [8]):  If  ip{t),  ao,  are  such  that 


+  00 


inf  r(®o7) 

1<  7  <ao  1 

-  -  j  =  -oo 


+  0O 


sup  E  mail) 


l<l7l<«0j=_oo 


>0, 


<  oo, 


(140) 


and  if 


OO 

Pijs)  =  sup  V  U(a^7) 

l<lTl<aOj=_oo 


^(ooT  +  7^) 


decays  as  fast  as  (1  +  with  e  >  0,  then  there  exist  6thr  >  0  such  that  the  ipj,k{t) 

constitute  a  frame  for  all  choices  bo  <  6thr-  For  bo  <  6thr,  the  following  expressions  are 
frame  bounds  for  the 


A 


inf 

l<|7|<ao 


E  Mr-EK^*)< 


k  =  —  oo 

k^O 


27/  . 


ll/2 
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{sup 

l<l7l<“0 


H-oo 


5^  +  XI 


2  r 


k—  —  oo 


The  condition  on  /?(7s)  and  Equation  140  are  satisfied  if  |■0(7)|  <  C'|7|“(l  +  |7|)  ^  with 
a  >  0,  id  >  a  +  1.  It  can  also  be  verified  that  if  the  '0j,fe(^)  Equation  131  constitute  a 
frame  with  frame  bounds  A  and  B  for  some  choice  of  Oo,  bo,  then  the  Fourier  transform  of 
ip{t)  satisfies 


boA  <  X  |^(‘^o7) 

j--CO 


<  boB 


(141) 


almost  everywhere  [74]. 

As  mentioned  before  it  is  advantageous  to  set  the  elementary  dilation  step  Oo  2 
because  it  leads  to  an  octave  by  octave  partitioning  of  the  frequency  domain.  Hereafter  we 
will  assume  this  partitioning  of  the  frequency  domain  unless  stated  otherwise. 


Multivoice  Wavelet  Frames 

In  Section  2.5.3  we  mentioned  that  it  is  advantageous  to  use  wavelet  frames  with 
H/A  -  1  <C  1,  because  the  can  be  approximated  by  'ipj,k{t)-  This  implies  that  the 

sum  in  Equation  141  should  be  almost  constant  for  7^0.  This  imposes  a  strong 
restriction  on  ip{t)  that  is  not  generally  satisfied.  In  order  to  overcome  this  problem 
without  having  to  give  up  too  much  freedom  in  choosing  V’(t)  or  its  bandwidth,  one  can 
adopt  the  use  of  different  voices  per  octave  (i.e.,  suboctave  sampling)  [8].  This  can  be 
achieved  by  using  different  wavelets  4>i{t) , 'ip2{t) t  and  looking  at  the  frame 

{A  -  kbo)  ,j,ke'L,v  =  l,...,N.\n  [8],  Daubechies  gives  the  following 
expressions  for  the  frame  bounds  of  this  multivoice  frame 


A 


B 


(  N  +0O  N 

E  A(2’7)P-EE 


j=-oo 


k^O  v=l 


N  H-oo 


N 


,  ,  sup  X  X  I^^(2^7)I^  +  XX 

1^171^2  k-^0  v=l 


with 

-f-OO 

A(7s)  =  sup  X  l^t>(2-^7)ll^t>(2^7  +  75)l- 
l<l7l<2j=_oo 
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It  is  suggested  in  [8]  that  by  choosing  the  ^'1(7))  •  •  ■  i'0Af(7)  to  have  slightly  staggered 
frequency  localization  centers,  coupled  with  good  decay  at  00,  one  can  achieve 
B/A  —  1  <C  1.  One  choice  favored  by  several  authors  [78,  8,  79]  is  to  take  “fractionally” 
dilated  versions  of  a  single  wavelet 

In  this  case  which  can  be  made  almost 

constant  by  choosing  N  large  enough  [8]. 

2.5.4  The  A  Trous  and  Mallat’s  Algorithms 

Recall  from  Section  2.5.3  that  the  wavelet  transform  of  a  signal  s{t)  is  given  by 

1  ft  —  b\ 

(W,ps){a,b)  =  la|“2  J  s{t)'if  dt,  s{t)  G  a, 6  €  M,  {a  ^  0). 

The  wavelet  coefficients  (s(t),  in  Equation  132  for  uq  =  2,  60  =  1  can  be  obtained  by 

sampling  the  wavelet  transform  of  s{t)  with  a  =  2^,  j  G  Z  and  b  =  2^k,  j,  A:  €  Z.  We  will 
refer  to  these  coefficients  as  the  decimated  wavelet  “series”  coefficients  and  will  denote 
them  by 

{C^s)j,k  =  {W^s){2\2^k) 

=  ^  (142) 

where  -  k)  is  the  generator  of  a  family  of  decimated  wavelets.  We  will 

refer  to  the  operation  of  computing  the  above  wavelet  coefficients  as  the  decimated  wavelet 
transform. 

For  many  applications  the  above  sampling  is  inadequate.  In  particular,  a  finer  grid  is 
needed  when  computing  the  wavelet  transform  for  applications  that  require  “good” 
time-frequency  localization.  This  can  be  achieved  by  computing  an  undecimated  wavelet 
transform  and  the  addition  of  voices  [66].  The  undecimated  wavelet  transform  of  a  signal 
s{t)  is  obtained  by  sampling  the  wavelet  transform  of  s{t)  with  a  =  2^,  j  ^  Z  and 
b  —  k,  k  ^  Z.  We  will  refer  to  these  coefficients  as  the  undecimated  wavelet  “series” 
coefficients  and  will  denote  them  by 

=  ms){2^,k) 

=  ,  (143) 

where  V’fc(^)  =  {^)  generator  of  a  family  of  undecimated  wavelets. 
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In  what  follows  we  will  review  two  separately  motivated  implementations  of  the 
wavelet  transform,  the  algorithme  a  trous  [13]  and  Mallat’s  [41]  multiresolution 
decomposition.  These  algorithms  are  both  special  cases  of  a  single  filter  bank  structure,  the 
discrete  wavelet  transform  [66],  the  behavior  of  which  is  governed  by  one’s  choice  of  filters. 
The  a  trous  algorithm  was  originally  devised  as  a  computationally  efficient  implementation 
of  the  wavelet  transform,  however  it  is  more  properly  viewed  as  a  nonorthogonal 
mnltiresolution  decomposition  for  which  the  discrete  wavelet  transform  computes  a 
“sampled”  continuous  wavelet  transform  (wavelet  “series”  coefficients)  exactly  [66].  On  the 
other  hand,  Mallat’s  algorithm  was  devised  as  an  orthogonal  multiresolution  decomposition 
[41]  that  when  initialized  properly  computes  a  “sampled”  continuous  wavelet  transform 
[66].  For  a  more  detailed  and  rigorous  account  of  these  algorithms  see  [41,  13,  66]. 

The  following  additional  terminology  and  notation  will  be  used.  The  discretized 
wavelet  transform  (with  sampling  period  T  =  1)  is  given  by 

(u;^s)(o,&)  = 

The  discretized  decimated  wavelet  “series”  coefficients  are  given  by 
{c^s)j,k  =  {w^s){2^,  2^k)  = 

^  ^  n 

Similarly,  the  discretized  undecimated  wavelet  “series”  coefficients  are  given  by 

{r^s)j,k  =  {w^s){2^,  k)  =  ^J2 

Signals  and  filters  in  bold  face  will  be  treated  as  vectors.  The  ^th  element  of  vector  / 
is  denoted  by  [f]k  =  fk-  The  symbol  f  will  be  used  for  the  mirror  filter  [f]k  =  fl  =  7-k- 
The  decimator  operator  is  denoted  by  a  matrix 

Dk,m  =  S{2k-m) 

~  ^2k,m 

where  5k, m  is  the  Kronecker  delta  and  S{k)  =  (5fc,o-  The  dilation  operator  is  denoted  by  a 
matrix 

Uk,m  =  6{k-2m) 

—  5k,2m- 

Convolution  of  f  and  s  is  denoted  by  [/  *  Convolution  followed  by 

decimation  becomes  [D{f  *  s)]/b  =  72m^k,m[f  *  s]m  =  [/  *  5]2A)  =  Xlm  hk-mSm-  The 
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following  shorthand  will  be  used  for  the  convolution  of  f  and  s,  [/  *  s]  =  Fs,  where 
Fm,n  =  fm-n-  It  is  easy  to  verify  that  [-D(/  *  s)]  =  DFs.  The  symbol  f  will  also  be  used 
for  the  Hermitian  transpose  of  matrices  [F^]m,n  =  fn-m  —  fm-n-  For  discrete  signals  the 
translation  operator  is  denoted  as 

~  ^m—ni 


and  for  continuous  signals  as 


{Trs){t)  =  -  r). 

The  2:  transform  of  a  discrete  signal  is  given  by 

n 

Finally,  we  state  the  following  two  identities  [66],  which  will  be  used  later.  For  any  F  of 
the  form  F^^n  =  fm-n  we  have  that 

[{DFyU  =  [{DF)%k-2in  (147) 


Y,[{DFy]n,kz'^  =  z^'^Hfiz^^).  (148) 

k  r=0 

The  Decimated  A  Trous  Algorithm 

The  a  trous  algorithm  was  originally  devised  as  an  efficient  implementation  for 
computing  wavelet  series  coefficients.  To  achieve  this,  the  wavelet  series  coefficients  in 
Equations  142  and  143  are  approximated  by  their  discretized  counterparts  in  Equations 
145  and  146,  respectively.  As  a  starting  point  consider  implementing  Equation  145. 

Clearly,  as  j  increases  ip{t)  must  be  sampled  at  progressively  more  points,  creating  a  large 
computational  burden  [13,  66].  The  solution  proposed  by  [13]  is  to  approximate  nonintegral 
points  via  a  finite  filter  ff  As  an  example,  let  be  the  filter  ;^(0.5, 1.0, 0.5).  Then, 

= 

k 


approximates  a  sampling  of  Let  g  be  a  filter  defined  by  Qn  —  ipi-n).  Then  the 

above  interpolation  can  be  computed  by  first  dilating  and  then  convolving  the  result 


J  '0(|),  for  n  even, 

I  H^(V) +  ^(^))  ’  lorn  odd, 

HI). 
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with  a  filter  which  leaves  the  even  points  fixed  and  interpolates  to  get  the  odd  points 
[66].  The  condition  that  f  be  the  identity  on  even  points  is  sufficiently  important  that 
merits  the  following  definition:  A  lowpass  filter  f  is  said  to  be  an  a  trous  filter  if  it  satisfies 

f2k  =  5{k)/^/2.  (149) 


The  \/2  is  simply  a  convenient  means  of  including  the  normalization  factor  of  Equation  145 
in  the  filter. 

It  can  be  verified  [66]  that  the  interpolation  operation  can  be  written  as  follows 


l/'»(V9'))„  =  [FWs'], 

^  V  / n~2k9k 
k 

This  result  and  Equation  145  leads  to  the  following  approximation  for  the  discretized 
decimated  wavelet  series  coefficients  at  =  1, 


{CqpS^  1,/C 


''  n 

~  ^  ''j  ^  ^  f  n— 2k— 2m9m 

n  m 

—  ^  ^  9k— m  ^  ^  f2m-n^n 

m  n 

=  [g*  D{f  *  s)]fe 

=  [GDFs]k, 


where  [s]„  —  Sn  =  s{n).  Proceeding  inductively,  one  can  find  the  following  approximation 
for  the  discretized  decimated  wavelet  series  coefficients  for  all  j  [66] 

~  [G{DFys]k.  (150) 

Remarkably,  the  right-hand  side  in  the  above  equation  is  the  discrete  decimated  wavelet 
transform  of  s  at  level  j  [66],  which  we  denote  as  follows^ 

[c,(s)1j  =  |G(r>F)'s]t.  (151) 

®The  notation  [cj\k  =  Cj^k  will  be  used  as  a  shorthand  for  [cj(s)]*. 
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Equations  150  and  151  yield  the  following  formula  for  the  discretized  decimated  wavelet 
series  coefficients 


{Ci,s)j,k  ^  [Cj(s)]fe.  (152) 

The  discrete  decimated  wavelet  transform  can  be  computed  iteratively  over  j  as  follows 

n 

=  [£>(/* 

Wj]k  —  ^  9k-n[^j]n 

=  [£?  *  (153) 

where  [so]k  =  Sfe-  Except  for  decimation  of  the  output  the  a  trous  algorithm  described  in 
[13]  is  given  by  Equations  152  and  153. 

The  original  a  trous  algorithm  made  no  statements  regarding  the  accuracy  of  the 
approximation  in  Equation  150  or  even  the  discretization  from  Equation  119  to  Equation 
144.  A  major  step  forward  towards  treating  this  question  lies  in  the  results  of  [66]  outlined 
in  Section  2.5.5. 

Mallat’s  Algorithm 

In  this  section  we  review  Mallat’s  multiresolution  algorithm.  For  a  more  detailed  and 
rigorous  account  of  this  algorithm  see  [41,  80,  66].  Mallat’s  algorithm  has  basically  the 
same  structure  as  Equations  153,  that  is, 

[sj+i]fc  =  [D{h*Sj)]k 

[dj+i]k  =  [D{g*Sj)]k  (154) 

where  [so]k  =  Sfc-  In  keeping  with  the  literature,  the  lowpass  filter  is  denoted  by  h  instead 
of  /.  The  constraints  on  h  and  g  which  ensure  an  orthonormal  multiresolution  analysis  are 
[41,  66]: 


1.  Perfect  reconstruction 


2.  Orthogonality  of  h  and  g  with  respect  to  even  shifts 

h2m-j92n-j  ~  0)  (1^^) 

j 

3.  Bandpass  condition  on  g 

=  (157) 

n 

4.  Lowpass  condition  on  h 

X^/!„  =  V2.  (158) 

n 

In  matrix  notation  Equations  155  and  156  may  be  written  as  follows 

{H^U)  {DH)  +  (G^U)  {DG)  =  I,  (159) 

{DH){G^U)  =  0,  (160) 

where  Hm,n  =  hm-n-  Multiplying  Equation  159  on  the  left  by  DH  and  using  Equation  160 
we  see  that 

{DH){H^U)=I. 

It  follows  that  h  and  its  shifted  versions  by  even  shifts  form  an  orthonormal  set 

3 

Similarly,  multiplying  Equation  159  on  the  left  by  DG  and  using  Equation  160  we  have 
that 

{DG){G^U)  =  I. 

It  follows  that  g  and  its  shifted  versions  by  even  shifts  also  form  an  orthonormal  set 

^  ^  92m-j92n—i  ~  ^m,nj  Ul,  U  G  Z. 
j 

These  two  results  together  with  Equation  156  imply  that  Equations  154  is  an  orthogonal 
decomposition  of  the  discrete  signal  Sj.  Furthermore,  Equation  155  implies  that 


From  Equations  155-158  it  follows  that  Equations  154  represents  a  wavelet 
decomposition  as  described  below  [41,  66]. 

Define  a  scaling  function  with  Fourier  transform  given  by 

OO  ^  _ 

r=l  ^ 

It  follows  that 

which  in  the  time  domain  takes  the  form 


^  -  k). 

k 

Therefore,  the  dilates  and  translates  of  ^(t), 

't'lAt)  =  ") 

have  the  property  [66] 

4^j  +  l,k{t)  —  ^  ^  ra0j,n(^)- 

n 

Define  a  wavelet  function  by 

i){t)  =  ^  -  k).  (161) 

k 

Then,  using  Equation  156  and  the  above  properties,  one  can  show  that  the  family  of 
wavelets 

are  orthonormal,  and  that  the  dj  are  the  coefficients  of  the  expansion  of  s{t)  in  terms  of 
'4’j,k{t)i  th^it  is 

(I -*=)*’ 

provided  that 

/+00  _ 

s{t)<j){t  -  k)dt. 

-OO 

Therefore,  Mallat’s  algorithm  when  initialized  properly  (see  the  equation  above)  computes 
exactly  a  “sampled”  continuous  wavelet  transform  [66]. 
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The  Undecimated  A  Trous  Algorithm 

As  mentioned  earlier  for  many  applications  the  time-frequency  sampling  grid  provided 
by  the  decimated  wavelet  transform  is  inadequate  when  “good”  time-frequency  localization 
is  required.  One  way  to  overcome  this  limitation  is  by  computing  an  undecimated  wavelet 
transform.  Following  the  same  approach  described  above,  the  wavelet  series  coefficients  in 
Equation  143  can  be  approximated  by  their  discretized  counterpart  in  Equation  146.  From 
Equations  145  and  146  it  follows  that  (c^s)yfe  =  (w^s){2^ ,2^k)  and  (r^s)y,fe  =  {w^s){2^,k) 
should  coincide  at  /t  =  0.  Utilizing  this  fact,  one  can  obtain  the  kth  undecimated  wavelet 
series  coefficient  (r^s)yfc  by  translating  the  signal  back  by  k  and  computing  the  discretized 
decimated  wavelet  series  coefficient  (c^(T_fes))j  g.  This  result  together  with  Equation  152 
yields  the  following  approximation  for  computing  the  discretized  undecimated  wavelet 
series  coefficients 


«  [cj{T_ks)]o,  (162) 

where  [s]„  —  Sn  —  s{n).  Remarkably,  the  right-hand  side  in  the  above  equation  is  the 
discrete  undecimated  wavelet  transform  of  s  at  level  j  [66],  which  we  denote  as  follows® 

[rj(s)]jt  =  [cj{T_ks)]o-  (163) 

Equations  162  and  163  yield  the  following  formula  for  the  discretized  undecimated  wavelet 
series  coefficients 

(rv,s)yfc  [r’i(s)]/fc-  (164) 

A  few  observations  about  the  decimated  and  undecimated  discrete  wavelet  transforms 
are  in  order.  Notice  that  in  general  Cj  is  not  translation  invariant 

n 

=  J2[G{DFy]k,n+m[s]n 

n 

^  Y^[G{DFy]k-m,n[s]n- 

n 

However,  if  one  replaces  m  by  2^m  in  the  above  equation  and  uses  Equation  147 ,  the  last 
step  becomes  an  equality  [66] 

[ffi'(^2Jjn^)]fc  =  [ffi(*)]fc-m-  (165) 

®The  notation  [rj]k  =  Vj^k  will  be  used  as  a  shorthand  for  [r_,  (s)]fc 
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Thus,  translating  s  by  2^m,  translates  octave  j  by  m. 

On  the  other  hand,  it  can  be  verified  that  Vj  is  translation  invariant 

['^  ~  [Cj  (Tm— /c^)]o 

Also,  from  Equations  163  and  165  it  follows  that  sampling  every  2^  points  produces 
exactly  cj^k,  that  is. 


Cjj^  —  Tj^23k’ 

An  iterative  algorithm  for  computing  the  discrete  undecimated  wavelet  transform  can  be 
obtained  by  taking  2:  transforms  [66].  From  Equations  151  and  163 

k  m  k 

=  Y^[G{DFy]„,,„z'“s(z). 

ra 

Applying  Equation  148  to  the  above  equation  it  follows  that 

k  r=0 

where  j  =  0  is  understood  to  mean  there  are  no  factors  of  /. 

It  is  easy  to  see  that  W  f  is  /  with  2^  -  1  zeros  inserted  between  every  pair  of  filter 
coefficients  and  that  its  2:  transform  is  f{z^').  Therefore  the  discrete  undecimated  wavelet 
transform  can  be  computed  iteratively  over  j  as  follows 

=  [iU^f)*sj)]k 
[rj]k  =  [{U^g)  *  Sj]k 

where  [sojfc  =  ^k-  The  above  result  together  with  Equation  164  is  essentially  the  a  trous 
algorithm  described  in  [13]. 

2.5.5  The  A  Trous  Algorithm  as  an  Exact  Wavelet  Transform 

In  the  previous  section  we  analyzed  two  separately  motivated  implementations  of  the 
wavelet  transform.  We  learned  that  the  a  trous  algorithm  was  originally  devised  as  a 
computationally  efficient  approximation  of  the  “sampled”  wavelet  transform  whose 


no 


implementation  can  be  viewed  as  a  discrete  wavelet  transform  for  which  the  filter  f  is  an 
interpolator  satisfying  Equation  149  and  the  filter  g  is  obtained  by  sampling  the  wavelet 
function  xp{t).  We  also  learned  that  Mallat’s  algorithm  is  a  discrete  wavelet  transform  with 
filters  h  and  g  subject  to  the  constraints  in  Equations  155-158  that  when  initialized 
properly  computes  exactly  a  “sampled”  continuous  wavelet  transform  for  the  wavelet 
function  'ipit)  in  Equation  161.  Next  we  will  review  the  results  of  [66]  which  show  that  the  a 
trous  algorithm  can  be  viewed  as  a  nonorthogonal  multiresolution  decomposition  for  which 
the  discrete  wavelet  transform  computes  a  “sampled”  continuous  wavelet  transform  exactly. 

Define  a  scaling  function  (^(t)  with  Fourier  transform  given  by 

00  ^  _ 


or  equivalently 


where 


(f>{t)  -  j^'^[{DF)\kV^x{2^  -  k), 


x{t)  = 


1  fort  e  [—1/2, 1/2), 
0  otherwise. 


I  UJ=0 


For  normalization  purposes  ^  =  \/2,  which  implies 

Y.h  =  V2 


(167) 


and 

/+00 

(j){t)dt  =  ^(0)  =  1. 

-OO 

Suitable  regularity  conditions  for  the  inverse  Fourier  transform  of  the  product  in  Equation 
166  to  converge  to  a  reasonable  behaved  function  may  be  found  in  [81,  82,  66].  One  can 
verify  [66]  that  the  dilates  and  translates  of  ^{t)  have  the  property 

[DF\k^n4^j^n  jt)  > 

n 

n 

Furthermore,  if  /  is  a  finite  filter,  then  0(t)  has  finite  support  [66]. 

Define  the  wavelet  function 


(^)  ~  E/  +  k)  —  ^  gl4){t  -  k). 

k  k 


(169) 
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The  recursion  in  Equation  168  implies  that 


Mi)  = 

n 

This  expression  along  with  Equation  169  yields 

4j(f)  =  ^[G(5p)ili,„^.o,„(t). 

n 

Using  the  above  expression  in  Equation  145  the  discretized  decimated  wavelet  “series” 
coefficients  take  the  form 

I 

I  n 

n  I 

=  J2[G{DFy]k,nY.s{im-n). 

n  I 

Define 

is)„  =  5;;;s(o^<r^,  (go) 

i 

then  the  above  equation  becomes 


which  is  the  discrete  decimated  wavelet  transform  of  the  discrete  signal  s„  defined  in 
Equation  170.  Therefore,  using  the  initialization  in  Equation  170  the  discrete  decimated 
wavelet  transform  computes  exactly  the  discretized  decimated  wavelet  “series”  coefficients 
{c^'s)j^k-  Similarly,  using  the  above  expression  for  Equation  142  we  see  that  the 

decimated  wavelet  “series”  coefficients  take  the  form 


i^ip'  ^)j,k  — 


/+00  _ 

■OO 

f"  s(t)  ( 5;;;iG{Df  )']t,  n^0,n(^)  I 

7-00  \  ^  / 

^[G(GF)i]>,„ 


'»+oo  _ 

:,n  / 


—  OO 
+  00 


/-t-OO  _ 

-  n)dt 

-OO 
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Using 


/4-00  _ 

s{t)4){t  -  n)dt  (171) 

■OO 

the  above  equation  becomes 

=  [G{DFys],, 

which  is  the  discrete  decimated  wavelet  transform  of  the  discrete  signal  defined  in 
Equation  171.  Therefore,  using  the  initialization  in  Equation  171  the  discrete  decimated 
wavelet  transform  computes  exactly  the  decimated  wavelet  “series”  coefficients 
Notice  that  170  is  the  discretized  version  of  Equation  171. 

Finally,  we  review  the  significance  of  the  a  trous  condition  in  Equation  149.  One  may 
show  that  a  filter  /  is  an  a  trous  filter  =  5{n)  [66].  Using  this  observation  it  can  be 

shown  that  if  /  is  a  trous  then  (1)  Equation  170  becomes  [s]„  =  s{n)  and  the  discrete 
decimated  wavelet  transform  of  the  sampled  signal  s(n)  yields  the  exact  discretized 
decimated  wavelet  “series”  coefficients  ic^'s)j,k  and  (2)  the  wavelet  "0  (t)  defined  in 
Equation  169  satisfies  i)' (n)  =  gl  [66]. 

In  summary,  given  f  satisfying  Equations  149  and  167,  g  defined  by  —  'ip{n)  (where 
7/;(t))  is  an  arbitrary  wavelet),  and  given  the  initialization  in  Equation  171,  the  discrete 
wavelet  transform  computes  an  exact  “sampled”  continuous  wavelet  transform  of  s{t)  using 
the  wavelet  function  ip'  {t)  defined  in  Equation  169.  Furthermore,  if  there  is  sufficient 
regularity,  ip' {t)  and  ip{t)  will  be  “close”  since  they  coincide  on  the  integers  up  to  the 
length  of  g. 

2.5.6  The  Sine-Gabor  Wavelet 

It  is  well  known  that  the  modulated  Gaussian  or  Gabor  function  is  the  only  function  that 
achieves  optimum  time-frequency  localization.  The  question  arises  as  to  whether  the  Gabor 
function, 

can  be  used  in  the  context  of  wavelet  analysis  and  synthesis.  Unfortunately,  this  is  not 
possible  because  the  Gabor  function  is  not  an  admissible  wavelet  (for  details  on  the 
admissibility  condition  for  wavelets  see  Section  2.5.3.)  Therefore,  optimum  time-frequency 
localization  via  the  wavelet  transform  is  not  possible.  However,  a  modified  version  of  the 
Gabor  function,  know  as  the  Morlet  wavelet  [8], 

ipM{t)  =  ^ 
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satisfies  the  wavelet  admissibility  condition,  and  for  cto  =  1,  wq  =  tt  (2/ In  achieves 
nearly  optimum  time-frequency  localization  and  generates  a  frame  of  wavelets  (for  more 
details  about  wavelet  frames  see  Section  2.5.3.)  In  fact,  the  error  incurred  in  approximating 
the  Morlet  wavelet  with  a  Gabor  function  under  these  conditions  is  negligible  [8].  Although 
Gabor  functions  are  not  wavelets  in  the  strict  sense,  the  term  “Gabor  wavelet”  is  used 
throughout  the  literature.  We  will  refer  to  Gabor  functions  as  Gabor  wavelets  provided 
that  is  negligible. 

Having  a  frame  of  wavelets,  j,  A;  e  Z,  is  a  desirable  property  because  it  allows 

one  to  completely  characterize  and  reconstruct  a  signal  s{t)  from  its  discrete  wavelet 
coefficients  (s(t), '0j,A;(A))  (see  Section  2.5.3.)  In  practice  it  is  desirable  to  have  a  frame  of 
wavelets  for  which  both  analysis  and  synthesis  can  be  computed  efficiently. 

We  show  that  the  sine-Gabor  function  (complex  part  of  the  Morlet  wavelet)  achieves 
nearly  optimum  time-frequency  localization  and  at  the  same  time  generates  a  frame  of 
wavelets.  Furthermore,  we  show  that  both  analysis  and  synthesis  using  the  sine-Gabor 
wavelet  can  be  computed  efficiently. 

The  fact  that  the  linear  combination  of  two  wavelets  is  also  a  wavelet  gives  us  an  easy 
way  to  show  that  the  sine-Gabor  function  is  a  wavelet.  By  taking  the  real  and  imaginary 
parts  of  the  Morlet  wavelet  we  see  that 

Retpuit)  =  (cos{coot)  -  , 

and 

Im  ipM{t)  =  sin(a;ot). 

It  follows  that  the  sine-Gabor  function  is  indeed  a  wavelet.  It  also  follows  that  the 
cosine-Gabor  function  is  not  a  wavelet  for  otherwise  the  Gabor  function  would  be  a  wavelet. 

The  question  as  to  whether  the  sine-Gabor  wavelet  generates  a  frame  of  wavelets  and 
the  efficient  computation  of  both  analysis  and  synthesis  using  the  sine-Gabor  wavelet  are 
addressed  below. 

2.5.7  The  Sine-Gabor  Wavelet  Frame 

We  showed  in  Section  2.5.6  that  sine-Gabor  functions  are  wavelets.  In  this  section  we 
study  the  time-frequency  localization  characteristics  of  the  sine-Gabor  wavelets  and  also 
establish  that  under  certain  conditions  they  generate  wavelet  frames. 

The  sine-Gabor  wavelet  is  given  by 

V>(t)  =  i^^e-"'/2-^sin(a;ot),  (172) 
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where 


K,!,  =  tt' 


-1/4 


(Jo(l  - 


1/2 


was  chosen  so  that  ||i/)||2  =  1.  The  time  and  frequency  resolution  of  the  sine-Gabor  wavelet 
are  given  by  the  root  mean  square  extent  of  the  function  and  its  Fourier  transform, 
respectively  (see  Equation  111  in  Section  2.5.2.)  It  follows  from  Equations  172  and  111 
that  the  time-resolution  of  the  sine-Gabor  wavelet  is  given  by 


o-(^) 


^o‘^o(l -2aga;g)) 
2(1  —  e~°'o“o) 


1/2 


The  Fourier  transform  of  the  sine-Gabor  wavelet  is  given  by 

g-<7g(a;-a)o)^/2  _  g-<7g(w-l-wo)^/2 
- 


where 


\/27r(Jo. 

Computing  the  root  mean  square  extent  of  '^(a;)  as  defined  in  Equation  111  would  lead  to 
misleading  results.  It  can  be  verified  that  if  a  function  '0(t)  satisfies  the  wavelet 
admissibility  condition  and  is  also  a  window  function  then  its  Fourier  transform  ^^(la)  must 
vanish  at  the  origin  (see  Equation  125  in  Section  2.5.3.)  In  addition  if  the  wavelet  i/>(t)  is 
real  we  have  that  ^(w)  =  ip{-oo),  so  that  |i/’(a;)|  is  a  symmetric  function  around  the  origin. 
It  is  therefore  more  suitable  to  restrict  our  attention  to  positive  frequencies  for  this  class  of 
wavelets  and  use 

[Mil  Jo 


for  the  frequency-resolution,  where 


It  is  easy  to  verify  that  the  sine-Gabor  wavelet  satisfies  all  the  conditions  mentioned  above 
and  that  its  frequency-resolution  is  given  by 


( 1  +  2olul  - 

2crg(l -e-<^o‘^o) 


1/2 
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Figure  35:  Square  of  the  product  of  the  time  and  frequency  resolution  for  the  sine-Gabor 
wavelet  as  a  function  of  (Tq  and  No¬ 


where 


m{ip)  — 


(Jo  erf((Toa;o) 


Figure  35  shows  the  square  of  the  product  of  the  time  and  frequency  resolution  for  the 
sine-Gabor  wavelet  as  a  function  of  ao  and  uq.  For  ctoo^o  >5/2,  the  sine-Gabor  wavelet  is 
within  1%  of  the  optimum  time-frequency  localization  lower  bound  given  by  the  uncertainty 
principle,  cr^('0)  >  1/4  (for  more  details  see  Equation  118  in  Section  2.5.2.) 

Next,  we  establish  that  the  sine-Gabor  function  does  indeed  generate  a  frame  of 
wavelets.  General  conditions  on  a  wavelet  -ipit)  under  which  a  frame  is  obtained  and 
estimates  for  the  frame  bounds  are  discussed  in  Section  2.5.3.  A  wavelet  'tp{t)  generates  a 
frame  of  wavelets  if  the  following  condition  is  satisfied 

|^(a;)|<C|n;r(l  +  la;|)-^  (173) 


with  a  >  0,  0  >  a  +  1- 

The  magnitude  of  the  Fourier  transform  of  the  sine-Gabor  wavelet  may  be  written  as 
follows 

|^(^)|  =  ^  (l  -  g-<72(|u;|-a;o)V2 

<  ^  (2(72a;o|w|)'/2e-<^§(l-'l-'o)V2 

^  (2aga;o|a;|)^/^ 

“  2  1 -h  cr^dwl  -  a;o)2/2’ 

where  the  first  and  second  inequalities  are  established  by  using 

1  -  Vo;  >  0,  0  <  7  <  1,  and  respectively.  It  follows  that  the 

sine-Gabor  wavelet  satisfies  Equation  173  with  a  =  0.4  and  (3  =  1.85  for  some  C 
sufficiently  large. 
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(c) 


Figure  36:  Sine-Gabor  wavelet  (left)  and  magnitude  of  its  Fourier  transform  (right)  with  (a) 
(70  =  5/2,  a;o  =  1,  (b)  (Jq  =  1,  wq  =  1  and  (c)  do  =  1.0657,  loq  =  0.0299. 


Frame  Bounds 

Table  6  and  Figure  36  show  the  estimated  frame  bounds  and  sine-Gabor  wavelet  for 
various  values  of  do  and  wo-  The  results  in  Table  6(a)  show  that  it  is  possible  to  obtain  a 
frame  of  wavelets  by  using  the  sine-Gabor  wavelet  and  at  the  same  time  achieve  nearly 
optimum  time-frequency  localization,  d^('0)d^('0)  =  0.2525.  Table  6(b)  shows  that  by 
trading  off  time-frequency  localization,  d2(^)d2(^)  =  0.3297,  the  frame  bounds  can  be 
made  tighter.  As  mentioned  in  Section  2.5.3,  snug  frame  bounds  are  a  desirable  property 
because  the  dual  wavelet  can  be  approximated  by  The  error  between  the  original 

signal  and  its  reconstruction  using  the  above  approximation  is  given  by  ^||/||2,  where 
r  =  B/A-l  (see  Section  2.5.3  for  more  details.)  For  the  sine-Gabor  wavelet  with  do  =  1, 
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Table  6(a).  Frame  bounds 

with  (To  =  5/2,  Wo  =  1 


bo 

A 

B 

BjA 

1.00 

4.3873 

5.6329 

1.2839 

1.25 

3.5099 

4.5063 

1.2839 

1.50 

2.9216 

3.7586 

1.2865 

1.75 

2.4089 

3.3169 

1.3769 

2.00 

1.6016 

3.4085 

2.1281 

2.25 

0.3963 

4.0572 

10.2385 

Table  6(b). 

Frame  bounds 

with  do 

=  1,  Wo  — 

1 

bo 

A 

B 

BjA 

0.50 

8.5076 

8.6988 

1.0225 

0.75 

5.6715 

5.7995 

1.0226 

1.00 

4.1918 

4.4114 

1.0524 

1.25 

2.8464 

4.0362 

1.4180 

1.50 

1.2635 

4.4720 

3.5395 

Table  6(c). 

Frame  bounds 

with 

do  —  1.0657,  Wo  = 

0.0299 

^0 

A 

B 

B/A 

0.75 

7.2030 

7.3285 

1.0174 

1.00 

5.3999 

5.4988 

1.0183 

1.25 

4.2456 

4.4733 

1.0536 

1.50 

3.1604 

4.1054 

1.2990 

1.75 

1.9203 

4.3075 

2.2432 

2.00 

0.6524 

4.7969 

7.3523 

Table  6:  Frame  bounds  for  the  sine-Gabor  wavelet  with  (a)  ao  =  5/2,  wq  =  1,  (b)  ao  =  1, 
Wo  =  1,  and  (c)  do  =  1.0657,  wq  =  0.0299.  The  elementary  dilation  parameter  ao  =  2  in  all 
cases. 


Table  7.  Frame  bounds 


with  cTo 

=  1.0657. 

bo 

A 

B 

B/A 

0.75 

7.2043 

7.3298 

1.0174 

1.00 

5.4008 

5.4997 

1.0183 

1.25 

4.2465 

4.4739 

1.0536 

1.50 

3.1615 

4.1055 

1.2986 

1.75 

1.9218 

4.3071 

2.2412 

2.00 

0.6541 

4.7962 

7.3328 

Table  7:  Frame  bounds  for  the  first  derivative  of  a  Gaussian  with  ao  =  1.0657. 
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Figure  37:  First  derivative  of  a  Gaussian  (left)  and  the  magnitude  of  its  Fourier  transform 
(right)  with  gq  —  1.0657. 

ljq  =  1,  the  error  is  2.55%  for  oq  =  2,bo  =  1.  The  parameters  for  the  sine-Gabor  wavelet 
used  in  Table  6(c)  were  obtained  by  using  a  Simplex  search  method  to  find  the  minimum  of 
the  difference  between  B  and  .4  as  a  function  of  gq  and  coq.  The  error  in  this  case  is  0.90% 
for  Co  =  2,6o  =  1-  This  reduced  error  is  at  the  expense  of  time-frequency  localization, 
a^(V’)o-^(^)  =  0.3401.  An  important  observation  follows  from  this  last  example:  For 
Go  =  1.0657,  coo  =  0.0299  the  sine-Gabor  wavelet  is  nearly  equaF  (up  to  a  sign)  to  the  first 
derivative  of  a  Gaussian  with  gq  =  1.0657.  Therefore,  the  first  derivative  of  a  Gaussian 
(Canny’s  approximation  to  the  optimal  step  edge  detector  [19])  can  be  used  to  completely 
characterize  and  reconstruct  a  signal  from  its  discrete  wavelet  coefficients.  Figure  37  shows 
the  first  derivative  of  a  Gaussian  and  the  magnitude  of  its  Fourier  transform.  Table  7  shows 
the  frame  bounds  for  the  first  derivative  of  a  Gaussian  with  gq  —  1.0657.  The  advantage  of 
the  sine-Gabor  function  over  the  first  derivative  of  a  Gaussian  is  that  by  changing  the 
product  of  GoUJo,  one  is  able  to  trade  off  between  Canny’s  criteria  for  designing  optimal  step 
edge  detectors  [19],  that  is,  the  product  of  signal-to-noise  ratio  and  localization,  SA,  and 
the  multiple  response  constraint,  r.  This  is  a  property  of  Canny’s  optimal  step  edge 
detector  that  is  not  possible  to  achieve  with  the  first  derivative  of  a  Gaussian.  Table  8 
shows  the  frame  bounds  for  the  sine-Gabor  wavelet  for  0.75  <  gqujo  <  3.75. 

2.5.8  Summary 

In  this  chapter  we  have  presented  the  tools  that  will  allow  us  to  compute  a  multivoice 
undecimated  wavelet  transform  exactly  for  an  arbitrary  wavelet  such  as  the  sine-Gabor 
wavelet.  Furthermore,  we  have  shown  that  the  sine-Gabor  wavelet  is  indeed  a  wavelet 
frame.  However,  we  found  that  having  nearly  optimum  time-frequency  localization  and 
obtaining  a  tight  frame  at  the  same  time  is  not  possible.  Both  constraints  are  important 
because  they  allow  for  “good”  time-frequency  localization  and  efficient  reconstruction, 

is  assumed  that  both  functions  have  the  same  L^-norm. 
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Table  8.  Frame  bounds  for 
0.75  <  (7o  <  3.75,  loq  =  1. 


CTo 

A 

B 

BjA 

3.75 

2.5933 

6.8524 

2.6423 

3.50 

2.9659 

6.5120 

2.1956 

3.25 

3.3385 

6.2107 

1.8603 

3.00 

3.7031 

5.9576 

1.6088 

2.75 

4.0546 

5.7629 

1.4213 

2.50 

4.3873 

5.6329 

1.2839 

2.25 

4.6920 

5.5627 

1.1856 

2.00 

4.9418 

5.5248 

1.1180 

1.75 

5.0866 

5.4618 

1.0738 

1.50 

5.0580 

5.2931 

1.0465 

1.25 

4.7875 

4.9433 

1.0325 

1.00 

4.1918 

4.4114 

1.0524 

0.75 

2.8851 

4.0837 

1.4154 

Table  8:  Frame  bounds  for  the  sine-Gabor  wavelet  for  0.75  <  a^oJo  <  3.75.  The  elementary 
dilation  parameter  ao  =  2  and  the  elementary  translation  step  bo  =  1  in  all  cases. 


respectively.  We  have  also  learned  that  a  better  (tighter)  frame  can  be  obtained  by  adding 
voices.  This  will  not  only  ease  the  reconstruction  but  is  also  a  requirement  for  applications 
in  time-frequency  analysis.  Our  current  work  is  concerned  with  extending  the  sine-Gabor 
wavelet  to  the  multivoice  framework  and  with  the  computation  of  multivoice  undecimated 
wavelet  transforms.  Our  future  work  includes  wavelet  reconstruction  in  the 
multivoice/undecimated  framework  as  well  as  extending  our  results  to  the  analysis  of  two 
dimensional  data  including  digital  mammograms. 

2.6  Circular  Mass  Recognition  Based  on  the  Hough  Transform 

2.6.1  Introduction 

The  Hough  Transform  (HT)  is  a  standard  method  for  shape  recognition  in  digital  images 
[83,  84].  It  was  first  applied  to  the  recognition  of  straight  lines  [85,  86]  and  later  extended 
to  circles  [87],  ellipses  [88]  and  arbitrarily  shaped  objects  [89].  Its  advantages  include 
robustness  to  noise,  robustness  to  shape  distortions  and  to  occlusions/missing  parts  of  an 
object.  Its  main  disadvantage  is  the  fact  that  computational  and  storage  requirements  of 
the  algorithm  increase  as  a  power  of  the  dimensionality  of  the  curve.  This  means  that  for 
straight  lines  the  computational  complexity  and  storage  requirements  are  O(n^),  for  circles 
0{n^)  and  for  ellipses  0(n^). 

We  begin  with  a  review  of  the  problem  of  circular  mass  recognition  using  the  HT. 
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Figure  38:  Demonstration  of  the  fact  that  the  center  of  circle  c,  O,  belongs  to  the  line  that 
perpendicularly  bisects  the  segment  defined  by  points  A  and  B  (both  points  belong  to  the 
circumference  of  the  circle). 

Even  though  there  have  been  attempts  towards  the  recognition  of  circles  via  the  standard 
3D  HT  [90],  it  has  been  recognized  that  there  is  a  need  for  a  decomposition  of  the  search 
space,  to  simplify  the  problem  both  in  terms  of  computation  and  storage.  Previous 
algorithms  have  been  based  on  the  property  that  the  normal  to  a  point  of  the  circumference 
of  a  circle  passes  through  the  center  of  the  circle  [91].  This  approach  works  well  for  high 
signal  to  noise  ratios  and/or  simple  environments.  As  the  signal  to  noise  ratio  decreases, 
the  accuracy  of  the  gradient  estimation  decreases  [92,  93].  The  fact  that  gradient— based 
methods  are  heavily  dependent  upon  the  accuracy  of  the  gradient  estimation  explains  why 
they  are  not  robust  to  noise.  Another  disadvantage  is  that,  sometimes,  the  edge  detector  of 
choice  does  not  provide  gradient  information.  A  comparative  study  of  various  HT  based 
techniques  for  circle  recognition  has  been  performed  by  Yuen  et  al.  [94]. 

The  approach  taken  in  our  study  was  to  decompose  the  3D  search  into  a  2D  HT  and 
ID  radius  histograming.  For  the  first  part,  instead  of  relying  on  a  2D  gradient-based  HT 
we  developed  a  2D  bisection  based  HT.  The  property  we  exploit  is  that  the  line  that 
perpendicularly  bisects  any  chord  of  a  circle  passes  through  its  center  (see  Figure  38). 
Experimentation  with  synthetic  data  demonstrated  that  our  approach  is  more  robust  to 
noise  than  gradient-based  techniques.  The  price  we  pay  for  robustness  is  an  increase  in  the 
computational  requirements  due  to  the  labeling  of  the  connected  components  that  our 
method  requires. 

The  second  part  of  the  algorithm,  ID  radius  histograming,  is  used  to  validate  the 
existence  of  identified  circles  and  calculate  their  radius.  We  show  that  extracting 
information  from  the  radius  histogram  is  not  a  trivial  task  and  we  devise  a  filtering 
technique  that  solves  this  problem. 

Next,  we  discuss  the  steps  of  our  algorithm.  In  this  section  we  give  details  about  the 
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implementation  of  the  algorithm  and  analyze  the  effects  of  various  factors,  such  as 
digitization  of  the  image,  discretization  of  the  parameter  space,  noise,  etc.,  on  the  accuracy 
and  the  computational  efficiency  of  the  method.  In  Section  2.6.3,  we  provide  sample 
demonstrations  of  the  algorithm  with  synthetic  and  real  images.  In  the  first  part  of  this 
section,  we  use  synthetic  data  to  demonstrate  the  robustness  of  the  technique  and  compare 
it  with  the  gradient-based  2D  HT.  In  the  second  part  of  this  subsection,  we  present  results 
of  processing  a  digitized  phantom  and  real  mammograms  from  the  Mammographic  Image 
Analysis  Society  (MIAS):  MiniMammography  Database.  Finally,  in  Section  2.6.5,  we 
present  some  conclusions  and  summary. 


2.6.2  Algorithm 

A  2D  Hough  Transform  based  on  bisection 

After  edge  detection,  the  resulting  connected  components  are  labeled.  Each  connected 
component  is  expressed  as  a  chain  of  the  coordinates  of  its  points.  We  connect  pairs  of 
points,  of  the  same  component,  using  a  sliding  window  (see  Figure  39).  If  the  coordinates 
of  the  points  are  A^xa.Va)  and  5(xb,?/b),  the  equation  of  the  line  that  perpendicularly 
bisects  AB  is 


y  ^  +  Vb 

Ua  -  Vb  ‘^{va  -  Vb) 


(174) 


All  points  (members  of  the  parameter  space)  belonging  to  this  line  have  their  votes 
increased  by  one.  Highly  voted  points  provide  an  indication  of  the  existence  of  circles. 
These  points  are  the  centers  of  such  circles. 

Since  we  are  dealing  with  digital  circles  the  points  of  their  circumference  are  affected 
by  digitization  and,  therefore,  do  not  exactly  satisfy  the  standard  circle  equation: 


r"  =  (a:  -  xof  +  {y-  (175) 

where  r  is  the  radius  of  the  circle  and  {xq,  yo)  are  the  coordinates  of  the  center  of  the  circle. 

Furthermore,  there  is  a  need  for  a  discretization  of  the  parameter  space  for  two  reasons: 

1.  Computational  efficiency.  It  is  impossible  to  account  for  all  the  circles  that  may  exist 
in  the  image,  and 

2.  The  line  that  perpendicularly  bisects  a  chord  of  the  circle  is  highly  unlikely  to  pass 
through  its  center.  Due  to  the  digitization  of  the  image,  pixels  belonging  to  a  digital 
circle  do  not  exactly  satisfy  Equation  175  and,  therefore,  their  chords  do  not  coincide 
with  the  chords  of  a  true  circle. 
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Figure  39:  Graphic  illustrating  our  implementation  of  a  2D  HT  for  circle  center  localization. 
Points  Po  and  P3  of  the  connected  component  Pq,  Pi,  •••  ,  P9,  are  assumed  to  belong  to 
the  same  circle  and  all  points  that  belong  to  the  line  that  perpendicularly  bisects  the  line 
segment  they  define  receive  votes.  The  same  process  is  repeated  for  the  pairs  Pi  P4,  P2  5 
etc.  For  this  particular  case  the  length  of  the  window  is  equal  to  three.  If  one  chooses  to 
connect  points  P0-P5,  Pi-Pe,  etc.,  the  length  of  the  window  is  equal  to  five. 


The  discretization  of  the  parameter  space  makes  the  voting  scheme  robust  to  noise  and 
errors  induced  during  the  edge  detection  and  other  preprocessing  operations.  On  the  other 
hand,  it  introduces  some  error  in  the  determination  of  the  position  of  centers  of  circles.  The 
larger  the  size  of  the  members  of  the  parameter  space,  the  higher  the  uncertainty  of  the 
estimation  of  the  position  of  the  center  of  the  circle.  In  our  implementation  the  parameter 
space  is  congruent  with  the  image  space  (this  is  dictated  by  the  adopted  parameterization) 
and  the  size  of  a  member  of  the  parameter  space  is  the  same  as  the  size  of  a  pixel.  This 
choice  is  reasonable  but  not  necessary.  For  example,  one  may  choose  a  more  robust  scheme 
where  the  size  of  a  member  of  the  parameter  space  is  2h  by  2h  {h  is  the  size  of  the  pixel). 
Such  a  discretization  will  not  only  increase  the  uncertainty  of  the  estimation  of  the  position 
of  the  center  of  the  circle  but  will  also  increase  the  probability  of  getting  accidental  peaks. 
Its  advantage  is  that  there  is  a  smaller  chance  of  missing  real  centers. 

Simple  detection  of  centers  of  circles  is  not  enough.  The  reason  is  that: 

1.  Highly  voted  pixels  provide  only  an  indication  of  the  existence  of  a  circle.  There  is  a 
need  for  verification  of  this  hypothesis.  Highly  voted  local  maxima  may  be  formed 
accidentally,  and 

2.  The  need  for  a  determination  of  the  third  parameter  of  the  circle  (its  radius  r). 

Next  we  describe  an  efficient  method  towards  the  verification  of  the  existence  of  a 
circle,  and  the  extraction  of  its  radius  via  radius— histograming. 
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Figure  40:  The  radius  histogram  for  the  center  of  a  circle.  All  bins  get  a  small  number  of 
votes  except  for  the  bin  to  which  the  circle  belongs.  A  sharp  maximum  appears  for  the  bin 
that  contains  the  circle. 

Analysis  of  radius  histogram 

After  detecting  possible  centers,  we  can  use  the  histogram  of  the  distances  of  all 
feature  points  from  the  centers  to  verify  the  existence  of  circles  and  extract  their  radii  (see 
Figure  40).  In  the  ideal  case  (continuous  circle,  exactly  determined  center)  the  analysis  of 
the  radius  histogram  would  be  an  easy  goal.  This  is  the  case  because  circles  would  show  up 
as  sharp  local  maxima  in  a  noisy  background  in  the  radius  histogram.  The  reasons  why 
this  does  not  happen  in  a  discrete  space  are: 

1.  Digitization/discretization  errors.  Even  for  the  case  where  we  can  exactly  determine 
the  center  of  the  circle,  the  digitization  of  the  circle  along  with  the  digitization  of  the 
image  combined  with  the  discretization  of  the  histogram  make  almost  certain  that 
the  votes  of  the  pixels  will  be  spread  to  a  number  of  neighboring  bins  of  the 
histogram.  As  one  can  see  from  Figure  41,  this  effect  is  similar  to  the  spreading  of 
the  straight  line  standard  HT  [95].  Under  the  boundary  quantization  scheme  [96],  it 
can  be  shown  that  digital  circles  can  be  bounded  by  two  Euclidean  circles  with  the 
same  center  and  radii  that  differ  by  h  [97].  This  is  illustrated  in  Figure  41  where  a 
digital  circle  is  bounded  by  continuous  circles  Ci  and  C2-  If  Ar  is  the  size  of  the  bins 
of  the  histogram  the  maximum  spreading  for  a  circle  is  equal  to: 

n,=  L^J+2,  (176) 

where  the  symbol  [z\  denotes  the  largest  integer  less  than  z, 

2.  Distortion  of  the  shape  due  to  imperfections  in  the  image  formation  (i.e.  breast 
compression),  errors  during  the  edge  detection  stage,  and  imperfections  in  the 
boundary  of  the  mass. 
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Figure  41:  All  pixels  belonging  to  a  digital  circle  can  be  bounded  by  two  concentric  circles, 
Cl  and  C2,  whose  radii  differ  by  h.  For  the  digital  circle  of  this  figure,  all  pixels  are  spread 
between  two  neighboring  bins,  bi  and  &2- 

3.  Noisy  pixels  and  other  objects  (tissues)  that  appear  in  the  neighborhood  of  the  circle. 

4.  Pixels  missing  from  the  boundary.  This  can  happen  because  of  occlusions,  missing 
parts  of  the  object,  etc. 

5.  Errors  in  the  localization  of  the  center  during  the  previous  steps,  mainly  during  the 
2D  Hough  voting  step.  If  the  localization  of  the  center  of  the  circle  is  not  accurate, 
two  peaks  will  appear  in  the  radius  histogram  of  the  estimated  center  (Figure  42). 
This  has  been  discussed  by  Yuen  et  al.  [94].  For  a  digital  image  if  the  error  in  the 
estimation  is  small,  say  less  than  two  pixels,  there  appears  a  single  extended  peak. 
The  length  of  the  peak  is  3  to  4  bins.  If  the  error  in  the  position  estimation  of  the 
center  is  larger,  two  local  maxima  appear  in  the  histogram.  As  the  error  increases  so 
does  the  distance  between  the  two  peaks.  Figure  42  illustrates  these  ideas. 

6.  The  number  of  pixels  belonging  to  a  digital  circle  and  a  digital  ring  increases  almost 
linearly.  Kulpa  [97]  showed  that  as  r  ->  +oo 

P,(r)  =  4\/2r,  (177) 

where  Pc  is  the  number  of  points  of  a  digital  circle  with  radius  r.  He  also  provided 
experimental  evidence  that  good  approximation  of  the  number  of  pixels  belonging  to 
a  digital  ring  (the  digital  object  bounded  by  two  continuous  circles  whose  radii  differ 
by  h)  is  given  by 

Priv)  =  27rr.  (178) 

We  should  emphasize  here,  that  the  equivalent  of  a  bin  in  a  digital  image  is  a  digital 
ring  and  not  a  digital  circle. 
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Figure  42;  If  an  inaccurate  estimate  of  the  center  of  the  circle  is  provided  by  the  first  step, 
the  radius  histogram  will  give  two  peaks.  In  this  example,  the  estimated  center  is  point  P, 
while  the  true  center  is  point  O.  Obviously,  these  peaks  will  appear  at  bins  with  distance  ri 
and  r2  from  the  estimated  center. 


All  these  factors  complicate  the  assessment  of  the  radius  histogram.  Our  approach 
towards  the  extraction  of  information  from  the  radius  histogram  is  based  on  the  previously 
presented  arguments. 

The  filter  we  propose  is  given  by  the  following  equation 
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(179) 


The  details  of  the  derivation  of  the  coefficients  of  this  filter  are  presented  in  Appendix  A.3. 
The  divisor  on  the  right  hand  side  of  Equation  179  is  equal  to  the  number  of  pixels 
belonging  to  a  circle  of  radius  r  (see  Equation  177).  We  should  add  that  the  proposed 
filtering  technique  gives  an  unbiased  estimator  of  the  normalized,  unoccluded  part  of  the 
circumference  of  the  circle.  In  other  words,  if  no  circular  features  exist  in  the  image  and  the 
image  is  corrupted  with  uniform  noise,  the  members  of  the  histogram  will  have  zero  mean. 

The  resulting  filtered  histogram  is  searched  for  local  peaks  whose  values  exceed  a 
certain  threshold.  Due  to  the  normalization  of  the  filtered  histogram  the  threshold  does 
not  depend  on  the  radius  of  the  circle. 


2.6.3  Results 

In  this  section  we  present  sample  results  from  experiments  applying  the  algorithm  with  real 
and  synthetic  digital  images.  Our  purpose  was  to  study  the  variation  of  robustness  for  the 
bisection-based  HT  with  increasing  noise  and  compare  it  with  gradient-based  HTs.  As  a 
measure  of  robustness,  we  used  the  robustness  ratio  P,  defined  by 

P  =  (180) 

'^total 
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(a) 


(b) 


Figure  43;  Comparison  of  the  robustness  P  as  a  function  of  the  SNR  for  three  circles  of 
radii:  (a)  r  =  20;  (b)  r  =  30;  and  (c)  r  =  50  pixels.  As  the  step  becomes  larger  P  increases. 


where  Vpeak  was  the  number  of  votes  of  the  peak  cell  and  Vtotai  was  the  number  of  pixels  of 
the  circumference  of  the  circle.  Obviously,  due  to  the  digitization  of  the  image,  even  in  the 
ideal  case,  P  is  less  than  1. 


Length  of  window 

The  first  thing  studied  is  the  dependence  of  P  with  respect  to  the  size  of  the  window 
slid  along  the  boundary  of  an  object.  Figure  39  shows  examples  of  windows  of  different 
sizes  slid  along  a  digital  curve.  Figure  43  shows  plots  of  P  for  three  circles  of  radii  20,  30, 
and  50  pixels.  We  observed  that  as  the  radius  increases  so  does  P.  We  also  noticed  that  as 
the  size  of  the  window  increases,  P  also  increased.  Any  theoretical  analysis  that  would 
indicate  the  optimum  window  size  is  difficult  due  to  the  fact  that  errors  in  the  coordinates 
of  the  two  extreme  points  of  the  window  are  not  independent  of  each  other.  Such  a  study 
has  been  presented  by  Amir  [98].  The  assumptions  this  study  relied  upon,  were 
independent  errors  in  the  coordinates  of  the  two  extreme  points  of  the  window  and  an  a 
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priori  knowledge  of  the  radius  of  the  circle.  Our  effort  in  the  application  of  digital 
mammography  is  towards  the  recognition  of  circles  of  varying  radii.  The  length  of  the 
window  should  be  kept  as  low  as  possible.  The  reason  for  this  is  the  need  for  recognition  of 
occluded  or  distorted  circles.  For  these  two  cases,  the  larger  the  window  the  less  points 
that  contribute  to  the  creation  of  peaks  in  the  parameter  space.  For  the  purpose  of  this 
research  we  chose  a  window  of  length  equal  to  twenty  as  a  reasonable  compromise  between 
robustness  to  noise  and  robustness  to  missing  (occluded)  regions. 

2.6.4  Synthetic  Images 
Image  Formation 

The  purpose  of  using  synthetic  images  is  to  compare  the  bisection-based  HT  with 
gradient-based  HT  in  terms  of  robustness  and  to  extract  conclusions  about  the  accuracy  of 
the  technique  as  a  function  of  the  Signal  to  Noise  Ratio  (SNR).  Our  effort  was  to  replicate 
the  image  acquisition  model  presented  by  Lyvers  and  Mitchell  [93].  According  to  this 
model  the  grey  level  value  of  a  pixel  is  given  by  the  following  equation 

p{k+0.5)Ax  p{l+0.5)Ay 

I{kJ)=  /  /  f{x,y)dxdy,  (181) 

J{k-0.5)Ax  J{l-0.5)Ay 

where  f{x,y)  denotes  a  continuous  image  and  I{k,l)  denotes  the  digital  image.  An  image 
containing  a  circle,  modeled  with  this  scheme,  is  presented  in  Figure  44(a).  One  can  easily 
notice  the  smoothness  of  the  edges  of  the  circle.  The  image  shown  in  this  figure  is 
corrupted  with  Gaussian  noise  of  standard  deviation  1.  The  background  value  is  equal  to 
120  and  the  foreground  140.  The  SNR  is  calculated  by  the  following  equation  [93]; 

SNR  =  20logio-  (182) 

a 

where  c  is  the  contrast  and  a  is  the  standard  deviation  of  the  noise. 

Examples 

Figure  44(b)  shows  the  output  of  edge  detection  and  thresholding  of  the  image  of 
Figure  44(a).  The  edge  detection  method  of  choice  was  the  one  proposed  by  Canny  [19].  In 
Figure  44(c)  we  present  the  resulting  voting  space  in  a  mesh  format.  The  coordinates  of 
the  identified  peak  were  (128, 128)  and  they  exactly  coincided  with  the  coordinates  of  the 
original  circle.  In  Figure  44(d)  we  show  the  radius  histogram  and  in  Figure  44(e)  we  show 
the  filtered  histogram.  Finally,  in  Figure  44(f),  we  superimpose  the  detected  circle  on  the 
original  image. 

The  same  process  was  followed  for  an  image  of  low  SNR  and  the  results  are  presented 
in  Figure  45.  We  note  here  that: 
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Figure  44:  Figure  illustrating  the  key  stages  of  the  algorithm  for  a  high  SNR  synthetic 
image  (SNR  ^26).  (a)  Original  image  with  a  circle  at  its  center.  The  radius  of  the  circle 
was  20  pixels;  (b)  edge  detected  image;  (c)  plotting  of  the  voting  space  in  a  mesh  format. 
The  peak  indicates  the  possible  existence  of  a  circle.  The  coordinates  of  the  peak  provide 
an  estimation  of  the  coordinates  of  the  center  of  the  circle;  (d)  histogram  of  the  feature 
pixels  as  a  function  of  their  distance  from  the  identified  center  (peak  in  the  voting  space); 
(e)  filtered  histogram.  The  identified  peak  was  at  r  =  20;  (f)  detected  circle  superimposed 
on  the  original  image. 
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1.  The  number  of  votes  of  the  center  is  smaller  in  the  low  SNR  case  than  in  the  high 
SNR  case,  and 

2.  The  spreading  of  the  votes  in  the  radius  histogram  is  higher  in  the  low  SNR  (more 
bins  share  the  same  votes) . 

To  counter  these  two  problems,  one  has  to  set  the  vote  threshold  of  the  2D  HT  to  lower 
values  and  to  use  a  larger  window  during  the  histogram  filtering  stage.  A  final  note  is  the 
fact  that  the  value  of  the  filtered  histogram  is  greater  than  one  is  attributed  to  the  votes 
due  to  noisy  pixels. 

Comparison  of  bisection-based  with  gradient-based  HT 

Robustness  in  the  context  of  noise,  complex  backgrounds  and  accuracy  in  the 
determination  of  the  center  of  the  circle  are  the  two  main  issues  of  concern  when  using  the 
bisection-based  2D  Hough  transform  for  center  detection.  We  use  the  robustness  ratio  F, 
defined  by  Equation  180,  and  the  distance  of  the  detected  center  and  the  true  center  as 
measures  of  the  performance  of  each  algorithm.  We  performed  experiments  for  a  wide 
range  of  SNRs  (2  to  26)  for  various  radii.  Experiments  were  repeated  for  20  times  and 
mean  values  are  plotted  in  the  figures  that  follow. 

In  Figure  46  we  present  plots  of  robustness  R,  as  defined  by  Equation  180  as  a 
function  of  the  SNR  for  circle  of  radii  20,  30,  and  50  pixels.  Comparing  the  R  of  the 
gradient-based  HT  and  the  bisection-based  HT,  we  can  conclude  that  the  former  is 
advantageous  for  high  SNRs,  while  the  latter  is  far  more  robust  for  low  SNRs. 

It  was  also  found  that  the  accuracy  in  the  determination  of  the  center  for  low  SNR 
with  the  bisection-based  method  did  not  decrease  considerably  while  the  gradient-based 
method  gave  large  errors  (see  Figure  47). 

RMI  phantom 

Figure  48(a)  shows  a  radiographic  image  of  an  RMI  156  (Gamex  Inc.,  Middleton,  WI) 
mammographic  accredited  phantom.  Three  types  of  mammographic  features  are  visible 
within  this  phantom  image,  simulating  difference  sized  circular  malignant  masses  as  well  as 
small  microcalcifications  and  fibers  of  varying  contrast.  Small  circumscribed  lesions  are  one 
of  the  more  important  signs  of  breast  cancer  which  may  be  detected  by  mammography  in 
asymptomatic  women.  The  detection  of  lesions  in  digital  mammograms,  is  one  of  the  main 
areas  digital  mammography  has  focused  on  [32,  99].  Mammographic  phantoms  like  the  one 
shown  in  Figure  48(a)  are  used  to  measure  the  performance  of  mammographic  systems  and 
a  mammography  unit  needs  to  be  able  to  visualize  at  least  four  fibers,  three  groups  of 
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Figure  45:  Figure  illustrating  the  key  stages  of  the  algorithm  for  a  low  SNR  synthetic  image 
{SNR  =  2).  (a)  Original  image  with  a  circle  at  its  center.  The  radius  of  the  circle  was  20 
pixels;  (b)  edge  detected  image;  (c)  plotting  of  the  voting  space  in  a  mesh  format.  The  peak 
indicates  the  possible  existence  of  a  circle.  The  coordinates  of  the  peak  provide  an  estimation 
of  the  coordinates  of  the  center  of  the  circle;  (d)  histogram  of  the  feature  pixels  as  a  function 
of  their  distance  from  the  identified  center  (peak  in  the  voting  space);  (e)  filtered  histogram. 
The  identified  peak  was  at  r  =  20;  (f)  detected  circle  superimposed  on  the  original  image. 
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(a) 


Figure  46;  Plots  of  the  mean  value  of  the  maximum  value  of  the  filtered  histogram  hmax, 
for  20  experiments,  as  a  function  of  the  SNR.  Continuous  lines  denote  results  using  the 
bisection-based  HT,  while  dashed  lines  denote  results  of  the  gradient-based  HT,  for  circles 
of  radii;  (a)  r  =  20;  (b)  r  =  30;  and  (c)  r  =  50.  While  for  high  SNR  (i.e.  SNR  >  20) 
the  gradient-based  HT  works  better,  for  low  SNR  drop  the  bisection-based  HT  gives  better 
results. 
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(a)  (b) 


Figure  47:  Plots  of  the  error  in  the  determination  of  the  center  of  the  circle  with  the  use 
of  the  bisection— based  method  (continuous  lines)  and  the  gradient— based  method  (dashed 
lines)  for  circles  of  radii:  (a)  r  =  20;  (b)  r  =  30;  and  (c)  r  =  50.  For  most  of  the  cases  the 
error  of  the  gradient— based  method  is  larger  than  the  error  of  the  bisection— based  method. 
As  shown  in  these  plots,  the  error  of  the  bisection-based  method  for  a  wide  range  of  noise 
levels  remained  smaller  than  a  pixel. 


133 


(c) 

Figure  48:  Figure  illustrating  the  steps  of  the  algorithm  on  a  digitized  phantom,  (a)  Image 
of  a  phantom;  (b)  output  of  Canny  edge  detection;  (c)  detected  circles  superimposed  on  the 
original  image. 

microcalcification  specks  and  three  masses  to  pass  the  American  College  of  Radiology 
accreditation  program  in  the  United  States.  The  size  of  the  image  is  512  by  512  pixels. 
Figure  48(b)  shows  the  results  of  edge  detection.  Eighteen  circle  centers  were  detected, 
when  the  voting  space  was  searched  for  local  maxima  that  had  more  than  30  votes.  After 
radius  histograming,  histogram  filtering  with  the  filter  Equation  179,  four  circles  were 
detected.  These  four  circles  are  superimposed  on  the  original  image  in  Figure  48(c). 

Mammogram 

Figure  49(a)  shows  a  512  by  512  image  that  was  cropped  from  a  digital  mammogram 
(mammogram  md005.pgm  from  the  MIAS  MiniMammography 
Database-http://slOd. smb.man.ac.uk/MIASmini.html).  The  image  contains  two 
overlapping  circular  masses  (biopsy  proven)  located  at  the  lower  middle  part  of  the  image. 
Figure  49(b)  shows  the  output  of  the  Canny  edge  detection,  while  Figure  49(c)  shows  the 
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circular  masses  detected  by  the  proposed  algorithm.  Both  masses  were  identified  but  as 
can  be  seen  from  the  same  figure  a  false  positive  also  appears  (circle  located  at  the  top 
middle  of  the  image). 

2.6.5  Summary 

We  presented  a  two-step  HT  for  the  detection  and  localization  of  circles  in  digital  images. 
The  first  used  a  bisection-based  2D  HT.  Experiments  with  synthetic  and  real  data  suggest 
that  the  bisection-based  HT  is  more  robust  to  noise  than  the  gradient-based  HT  while 
being  more  accurate.  Davies  was  the  first  to  use  this  property  for  the  localization  of  centers 
of  circles  [91].  His  implementation  is  very  efficient  in  terms  of  computation  but  since  it 
used  only  vertical  and  horizontal  chords,  it  lacked  robustness  needed  for  analysis  in 
complex  backgrounds. 

The  second  step  used  radius  histograming  to  detect  a  circle  and  extract  its  radius. 
Enhancement  of  local  maxima  with  the  use  of  a  Laplacian  filter  as  suggested  by 
Kierkegaard  [100],  had  a  significant  shortcoming  that  it  did  not  normalize  for  the 
increasing  number  of  pixels  belonging  to  a  bin  when  moving  away  from  a  center  of  the 
circle  (see  Equations  178  and  177).  Furthermore,  there  was  no  intuitive  method  for  the 
selection  of  the  threshold  of  the  filtered  histogram.  Our  filtering  scheme,  considers  votes 
due  to  noisy  pixels,  shape  distortions  and  normalizes  to  account  for  the  dependence  of  the 
number  of  pixels  on  the  radius  of  a  circle. 

Sample  results  were  provided  and  showed  that  the  method  is  capable  of  detecting 
circular  masses  of  varying  sizes  and  shape  distortions  in  complex  background,  including 
dense  breast. 
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Figure  49:  Figure  illustrating  the  steps  of  the  algorithm  on  a  digitized  mammogram,  (a) 
Digitized  mammogram;  (b)  output  of  Canny  edge  detection;  (c)  detected  circles  superim¬ 
posed  on  the  original  image.  The  algorithm  correctly  detected  the  two  overlapping  masses 
while  giving  a  false  positive  (top  middle  of  the  image). 


136 


3  Conclusions 


During  the  past  year,  we  have  made  significant  progress  in  the  development  of  a 
methodology  for  accomplishing  adaptive  contrast  enhancement  by  multiscale 
representations.  Our  studies  have  demonstrated  that  features  extracted  from 
multiresolution  representations  can  provide  a  mechanism  for  the  local  emphasis  of  salient 
and  subtle  features  of  importance  to  mammography.  The  improved  contrast  of 
mammographic  features  makes  these  techniques  appealing  for  computed  aided  diagnosis 
(CAD)  and  screening  mammography.  Screening  mammography  examinations  are  certain  to 
grow  substantially  in  the  next  few  years,  and  analytic  methods  that  can  assist  general 
radiologists  in  reading  mammograms  shall  be  of  great  importance. 

In  the  paragraphs  below,  we  summarize  our  progress  and  identify  future  directions  of 
research  to  be  carried  out  during  the  final  year  of  our  investigation. 

The  one-dimensional  discrete  dyadic  wavelet  transform  was  extended  to  higher-order 
derivatives  and  even-order  spline  functions,  and  an  improved  initialization  procedure  was 
developed.  In  comparison  to  the  originally  employed  initialization  [7],  our  method  showed 
significantly  better  performance  at  finer  scales  of  analysis  (of  importance  for 
mammographic  features  including  microcalcifications). 

A  multiscale  spline  derivative-based  transform  was  constructed  as  an  approximation 
to  a  steerable  dyadic  wavelet  transform.  The  transform  used  x-y  separable  filters  in  a 
perfect  reconstruction  filter  bank  and  enabled  fast  directional  analysis  of  images  without 
the  introduction  of  artifacts  due  to  translation  and  rotation  invariance.  Such  artifacts  are 
inherent  to  traditional  methods  of  wavelet  analysis. 

Preliminary  image  fusion  results  showed  promise  for  producing  enhanced  images  with 
a  more  familiar  appearance  to  radiologists. 

We  have  described  a  parallel  algorithm  to  accomplish  high-speed  interactive  multiscale 
processing  for  enhancement  within  ROIs  in  digital  mammograms.  The  implementation 
relied  upon  foundation-level  libraries,  XIL,  on  top  of  a  SPARCstation’s  SX  frame  buffer. 
We  showed  the  amount  of  speedup  for  the  method  compared  to  traditional  linear 
convolution  and  a  conventional  FFT  approach.  Our  multiscale  approach  employing 
splitting  in  decomposition  and  merging  in  reconstruction  was  shown  to  be  efficient. 

We  showed  that  the  continuous  wavelet  transform  of  a  band-limited  signal  can  be 
estimated  from  samples  of  a  signal.  If  a  mother  wavelet  is  compact  in  the  frequency 
domain,  effectively  no  approximation  error  was  observed  experimentally. 

We  described  a  modified  dyadic  wavelet  transform,  which  was  able  to  decompose  a 
signal  onto  an  arbitrary  scale  and  reconstruct  perfectly  the  signal  from  that  scale.  We 
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presented  a  special  class  of  spline  wavelet  filters,  which  allowed  us  to  use  a  high  order 
scaling  function  without  significantly  increasing  the  computational  complexity  of  analysis. 

The  importance  of  arbitrary  scale  analysis  was  demonstrated  by  digital  radiographs  of 
a  mammography  phantom  and  digitized  mammograms.  The  study  showed  that  the 
algorithm  was  able  to  detect  very  subtle  masses,  which  were  rated  to  be  almost  invisible  by 
radiologist  specializing  in  mammography. 

An  enhancement  algorithm  relying  on  multiscale  wavelet  analysis  and  extracting 
oriented  information  at  each  scale  of  analysis  was  investigated.  The  evolution  of  wavelet 
coefficients  across  scales  characterized  the  local  shape  of  irregular  structures.  Using 
oriented  information  to  detect  the  features  of  an  image  appears  to  be  a  promising  approach 
for  enhancing  complex  and  subtle  structures  of  the  breast.  Steerable  filters  which  can  be 
rotated  at  arbitrary  orientations  can  reliably  find  visual  cues  within  each  spatial-frequency 
sub-band  of  an  image.  “Coherence  measure”  and  “dominant  orientation”  clearly  helped  us 
to  discriminate  features  from  complex  surrounding  tissue  typical  in  mammograms. 

Existing  and  previous  multiscale  enhancement  approaches  [72,  32,  59]  attempted  to 
enhance  an  image  by  detecting  edges.  Unfortunately,  most  edge  detection  algorithms  can 
not  distinguish  between  “authentic”  edges  and  phantom  edges.  In  contrast,  this  algorithm 
relied  upon  a  coherence  map  and  phase  information  which  resulted  in  an  enhancement 
naturally  close  to  the  original  image.  This  type  of  artifact  free  enhanced  image  can  provide 
more  “obvious”  (and  familiar)  visnal  cues  for  radiologists. 

In  this  chapter  we  have  presented  the  tools  that  will  allow  us  to  compute  a  multivoice 
undecimated  wavelet  transform  exactly  for  an  arbitrary  wavelet  such  as  the  sine-Gabor 
wavelet.  Furthermore,  we  have  shown  that  the  sine-Gabor  wavelet  is  indeed  a  wavelet 
frame.  However,  we  found  that  having  nearly  optimum  time-frequency  localization  and 
obtaining  a  tight  frame  at  the  same  time  is  not  possible.  Both  constraints  are  important 
because  they  allow  for  “good”  time-frequency  localization  and  efficient  reconstruction, 
respectively.  We  have  also  learned  that  a  better  (tighter)  frame  can  be  obtained  by  adding 
voices.  This  will  not  only  ease  the  reconstruction  but  is  also  a  requirement  for  applications 
in  time-frequency  analysis.  Our  current  work  is  concerned  with  extending  the  sine-Gabor 
wavelet  to  the  multivoice  framework  and  with  the  computation  of  multivoice  undecimated 
wavelet  transforms.  Our  future  work  includes  wavelet  reconstruction  in  the 
multivoice/undecimated  framework  as  well  as  extending  our  results  to  the  analysis  of  two 
dimensional  data  such  as  images. 

Finally,  we  presented  a  two-step  HT  for  the  detection  and  localization  of  circles  in 
digital  images.  The  first  used  a  bisection-based  2D  HT.  Experiments  with  synthetic  and 
real  data  suggest  that  the  bisection-based  HT  is  more  robust  to  noise  than  the 
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gradient"based  HT  while  being  more  accurate.  Davies  was  the  first  to  use  this  property  for 
the  localization  of  centers  of  circles  [91].  His  implementation  is  very  efficient  in  terms  of 
computation  but  since  it  used  only  vertical  and  horizontal  chords,  it  lacked  robustness 
needed  for  analysis  in  complex  backgrounds. 

The  second  step  used  radius  histogramming  to  detect  a  circle  and  extract  its  radius. 
Enhancement  of  local  maxima  with  the  use  of  a  Laplacian  filter  as  suggested  by 
Kierkegaard  [100],  had  a  significant  shortcoming  that  it  did  not  normalize  for  the 
increasing  number  of  pixels  belonging  to  a  bin  when  moving  away  from  a  center  of  the 
circle  (see  Equations  178  and  177).  Furthermore,  there  was  no  intuitive  method  for  the 
selection  of  the  threshold  of  the  filtered  histogram.  Our  filtering  scheme,  considers  votes 
due  to  noisy  pixels,  shape  distortions  and  normalizes  to  account  for  the  dependence  of  the 
number  of  pixels  on  the  radius  of  a  circle. 

Sample  results  were  provided  and  showed  that  the  method  is  capable  of  detecting 
circular  masses  of  varying  sizes  and  shape  distortions  in  complex  background,  including 
dense  breast. 

These  results  are  encouraging  and  continue  to  support  that  wavelet  based  image 
processing  algorithms  can  play  an  important  role  in  improving  the  imaging  performance  of 
digital  mammography.  We  believe  we  are  now  close  to  our  original  schedule  and  anticipate 
that  we  shall  complete  the  objectives  described  in  Phase  IV  and  Phase  V  this  year. 
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A  Appendix 

A.l  Initial  Condition 


Figure  50:  Modify  initial  condition 


A  non-ideal  acquisition  device  is  shown  in  Fig  50  (a).  The  impulse  response  of  the 
acquisition  device  is  ^{t),  where  ^{t)  is  a  lowpass  filter.  The  approximation  of  ^{t)  in 


space  is 


^{t)  «  |(t)  =  z{k)(j){t  -  k), 


(183) 


where  z(k)  is  the  projection  of  ^{t)  onto  the  scaling  space.  The  discrete  signal  x{k)  if 

/OO 

f{t)^{n-t)dt 

-OO 

/OO 

k 

/OO 

f{t)^(n  -  t  -  k)dt) 

k 

where  x'{n)  =  f{t)(p{n  —  t)dt  .  Given  x(k), 

x'(n)  =  —  k). 


(184) 


(185) 


When  ^(t)  is  an  ideal  lowpass  filter,  it  can  be  approximated  as  a  cardinal  spline  filter  [45]. 
If  the  scaling  function  is  a  spline  of  order  n  -  1, 


«*) «  E  ('>”  )  ^{k)  beta{t  —  k). 


(186) 


Therefore,  z-^{k)  =  6"(A:)  ,  z{k)  =  (5”)  ^(/c).  Notice  that  Equation  (185)  is  similar  to 
projecting  a  signal  onto  a  scale  space  as  described  in  Unser’s  algorithm  [43]  [44],  although 
the  starting  point  is  quite  different  in  our  case. 
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A. 2  Proof  of  Theorem  3  in  Section  2.3.7 


In  this  appendix,  we  will  show  that  if  (^(t)  =  /?”  ^(t),  4>i{t)  =  /?”*  ^{t)  , 


Pa{oo) 


^LO) 


Rect{uj)  )  *  2-!r52TT{co) 


(187) 


as  n  +  m  ^  +oo. 

Proof:  Note  a  central  B  spline  is  symmetric,  In  the  oblique 

projection  method, 

/oo  roo 

-  k)dt  =  /  -  t)dt  =  *  r-" 

-oo  'J  — oo 

Substituting  Equation  (188)  into  Equation  (85), 

/OO 

-  t)dtw*<"-')-\k  -  i) 

■OO 


(188) 


iez 

r^OO 


/OO 

-  i)dt. 
iez 


Since  in  the  frequency  domain,  we  write 


(189) 


(190) 


;  n— 1 


Define  (5  (w)  =  Then 


-t)  =  {i  -  t) 

/OO 

^"■'■™'“i(r))/i"“i(i  —  t  —  v)dv 

-OO 


let  u  =  i  —  V 

r*00 


/OO 

pn+m-l^j^  _  _  t)du. 

■OO 

Substituting  Equation  (191)  into  Equation  (189), 

/OO  poo 

-  «)/?“■'{«  -  i)rf»)((>“+’”-')-'(*:  -  i)dt 

■°°  iez  d-co 

nco 

-  u)P^-\u  -  -  i)dtdu 

■°°  iez 
poo  poo 

=  /  (  /  ^a{tW-\u  -  t)dt){Y,ld"^'"~\i  -  uW^^-^)-\k  -  i))du 

J-co  J-OO  jg2: 

Let  i'  =  k  —  i,  and 

poo  poo 


(191) 


-OO  •/  — oo 
r»oo  /‘OO 


i'ez 


(192) 


-OO  O  — OO 
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Where  ’7"^"*  ^(^)  ^  cardinal  spline 

filter.  As  n  +  m  +oo  ,  it  converges  to  the  ideal  sine  interpolator.  Therefore, 


Pa{0j) 


Rect{(jj) 


*  27f62niu}). 


(193) 


A. 3  Derivation  of  Equation  179 

In  this  appendix  we  will  derive  the  filter  shown  in  Equation  179.  We  shall  assumpe  that  we 
are  dealing  with  an  image  corrupted  with  uniform  noise  and  that  we  are  looking  for  a 
perfect  digital  circle.  As  we  noted  (see  Equation  176)  earlier,  for  Ar  =  h  the  maximum 
number  of  bins  a  circle  spreads  its  votes  is  equal  to  2.  Therefore,  if  we  pass  the  radius 
histogram  with  the  filter 

t=[l  1] 

we  will  obtain  an  estimate  of  the  pixels  the  circle  contains.  We  would  like  the  filter  to  be 
symmetric  around  the  central  bin  and  therefore  we  choose 

t  =  [i  1  1]. 

The  resulting  array  contains  the  number  of  pixels  of  the  circle  plus  votes  due  to  noise, 
other  objects  etc.  We  need  to  enhance  our  scheme  with  a  module  that  estimates  the  votes 
due  to  noise  and  subtracts  them.  For  this  purpose  we  use  two  more  bins  (the  outer  ones). 
Therefore  our  filter  becomes 


t  =  [a  1  1  1  P].  (194) 

To  get  values  for  a  and  /?  we  will  require  that  our  filter  give  an  unbiased  estimator 
with  minimum  uncertainty  (standard  deviation).  We  assume  that  noise  is  spatially 
uniform.  To  have  an  unbiased  estimator  means  that  if  the  image  is  corrupted  with  noise 
and  no  circular  object  is  present,  the  resulting  (processed)  histogram  will  have  mean  value 
equal  to  zero. 

It  is  easy  to  see  that  the  r-th  bin  of  the  histogram  represents  a  ring  [97]  whose  radius 
is  equal  to  r.  Kulpa  provided  experimental  evidence  that  a  good  approximation  to  the 
number  of  pixels  belonging  to  such  a  ring  is  27rr.  This  estimation  is  based  on  empirical 
evidence  [97]  and  until  now  no  theoretical  justification  to  it  had  been  found.  We  verified, 
through  simulations,  that  this  formula  holds  with  adequate  accuracy  for  radii  in  the  range 
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[h,512h].  It  is  easy  to  show,  that  the  number  of  votes  the  r-th  bin  takes  due  to  noise 
follows  a  binomial  distribution  with  mean  value: 

yU  =  p(27rr)  (195) 

and  standard  deviation: 

a  =  \/p(l  -p)(27rr),  (196) 

where  p  is  the  probability  that  a  pixel  is  noisy. 

From  Equation  195  and  the  constraint  we  imposed,  that  our  estimator  should  be 
unbiased,  we  formulate 

Pf  —  -  2h)  +  (r  -  h)  +  (r)  +  {r  +  h)  +  P{r  +  2h)  —  0.  (197) 

The  second  equation  that  we  use  comes  from  the  requirement  of  minimization  of  the 
standard  deviation  of  the  resulting  filtered  histogram.  The  standard  deviation  of  the 
filtered  space  is  given  by  the  following  formula; 

ap  =  2^ y/a^r  -  2hY  +  (r  -  hf  +  {rf  +  (r  +  hf  +  /^^(r  +  2hf.  ( 198) 

If  we  make  the  following  substitutions; 

XI  =  a{r  —  2h) 

X2=r-h 
X3  =  r 
XA=r+h 
X5  =  f3{r  +  2h) 

to  the  Equations  197  and  198  we  obtain 

Xl  +  X2  +  X3  +  XA  +  X5=-0,  (199) 

and 

ap  =  27r\/xl2TX22TX32TX42T^.  (200) 

If  we  solve  Equation  199  for  XI,  and  substitute  into  Equation  200,  and  differentiate 
Equation  200  with  respect  to  (5  we  obtain 

^  =  (2X5  +  2(X2  +  X3  +  X4  +  X5))/  =  0,  (201) 

op 

where  /  is  a  factor  that  cannot  become  zero  (a  constant  divided  by  a  function  of  X5). 

From  Equations  197  and  201  we  can  compute  a  and  {3.  The  resulting  filter  was 
described  in  the  main  part  of  Section  2.6  (see  Equation  179). 
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